Inversion for sound speed profile in shallow water based on long short-term memory networks and ray theory
-
摘要:
针对水声多途信道环境下的水下声速剖面反演问题, 将深度学习与射线声学理论结合, 提出了一种基于长短期记忆网络(LSTM)的反演方法。该方法根据有序线阵等距的特点, 将到达时间差、到达角度等多模态数据融合组成的感知矩阵作为输入, 利用LSTM网络处理时序数据的能力挖掘空间上有序分布的接收阵元之间的关联信息, 从而实现声速剖面反演。在此基础上, 还提出了基于互相关函数的硬阈值估计方法, 通过降低感知矩阵的测量误差提高模型的抗多途性能。通过数值仿真验证了该反演方法的可行性和准确性, 与传统优化算法相比, 所提算法能够更好地捕捉声速剖面的非线性特征, 具有更高的反演精度和较强的抗噪能力。
Abstract:To address the problem of underwater sound speed profile (SSP) inversion in underwater acoustic multipath channels, this paper combines deep learning and ray theory to propose an inversion method using a long short-term memory network (LSTM). Based on the equidistant characteristics of the horizontal line array, the proposed method takes the perceptual matrix composed of multi-modal data, such as time difference of arrival and angle of arrival, as input, and utilizes the ability of LSTM network to process time-series data to mine the correlations between spatially ordered receiving array elements for sound speed profile inversion. On this basis, a time delay estimation method based on hard threshold estimation method and cross-correlation function is proposed to reduce the measurement errors of the perceptual matrix and improve the anti-multipath performance. The feasibility and accuracy of the proposed method are verified through numerical simulations. Compared with the traditional optimization algorithm, the proposed algorithm better captures the nonlinear characteristics of SSP, with higher inversion accuracy and stronger noise resistance.
-
引言
声速剖面(SSP)在海洋声学研究中扮演着至关重要的角色, 其具备的时空维度上的动态特性会显著影响声波在水下的传播轨迹, 进而改变传统探测模型结构, 引入算法层面的误差。由此可见, 声速剖面的正确与否直接影响水下目标探测结果的精度和可靠性, 因此, 声速剖面的获取对于优化水下声探测、水下声定位、声呐性能评估以及声呐系统设计等方面具有极其重要的意义。然而, 如何准确地估计声速剖面一直是一个极具挑战的研究课题。目前, 获得声速剖面的方法主要包括直接测量法, 基于经验公式的间接测量法以及特征匹配法。直接测量法通常使用声速剖面仪进行测量。然而, 在使用过程中, 由于传感器受环境因素和自身漂移的影响, 声速剖面仪的稳定性会下降, 测量精度会降低, 且使用成本高昂[1]。间接测量法根据海水中温度、盐度和深度等参数, 通过声速经验公式计算得到海水声速。然而, 不同经验公式有不同的使用范围, 因此该方法不适用于大范围水域的声速剖面测量[2]。除了上述方法外, 还可以通过特征匹配法获得声速剖面。特征匹配法中最经典的方法是匹配场处理(MFP), 其通过水下声传播物理模型生成的特征参数与实际测量的特征参数的匹配反演环境参数。文献[3]利用匹配波束构造了反演代价函数, 并在此基础上研究了浅海声速剖面的匹配波束反演方法。文献[4]采用MFP技术与经验正交函数(EOF)分解相结合的方法, 成功实现了声速剖面反演。相较于前两种方法, 特征匹配的方法虽然具备实现简单、可靠性高、灵活性强以及适用范围广等优点, 但也存在依赖模型、特征参数选择困难等问题。
特征匹配中常用的物理模型有射线声学模型和简正波模型[5], 常用于匹配的特征参数主要有信号飞行时间[6]、信号模态相位[7]、简正波水平折射角[8]、信号声压[9]以及信号波束[10]等。与其他信号特征相比, 使用飞行时间作为特征参数反演声速剖面具有诸多优势。首先, 它受海底声速、密度等参数影响较小, 能够解决传统方法中边界参数不匹配的问题; 其次, 与其他特征量相比, 本征声线的直达飞行时间较容易获得; 最后, 由于飞行时间较容易获取, 三维空间下的声速剖面反演可以通过测算分布式节点飞行时间的方式实现。文献[11]介绍了基于飞行时间的声层析方法的实现, 通过物理模型建立飞行时间与声速场的联系进而得到反演式。文献[12]将声传播模型从非线性优化问题转化成线性方程组的形式, 并结合遗传算法实现了基于飞行时间的声速剖面反演。文献[13]研究了利用信号在水下固定参考点与水面已知位置之间的飞行时间来反演海水声速剖面的方法, 利用牛顿迭代法求解逆问题, 建立一种针对目标海域进行快速、精确的声速剖面反演方法。
尽管使用飞行时间反演声速剖面具有诸多优势, 但这种方法存在一些不可避免的缺点。一方面, 水声通信的收发端存在同步问题, 这导致飞行时间的获取存在同步误差, 虽然通过双向通信能够一定程度上消除误差, 却带来双倍的功率开销。另一方面, 水声信道具有严重的路径衰减、噪声、有限带宽和严重多途等特点, 这导致获得的声信号产生畸变, 进而使得获取到的飞行时间存在测量误差。针对这个问题, 文献[14]提出使用到达时间差(TDOA)作为特征进行反演, 但到达时间差对时延估计方法的误差较为敏感, 该文献并没有对其展开详细讨论。
此外, 所用反演方法的鲁棒性和泛用性也存在问题。在海况复杂的海域, 测量误差较大的特征会导致反演模型的性能恶化。为了解决这些问题, 研究人员开始尝试使用机器学习方法进行水下声速剖面反演。这种方法具有从数据集中提取非线性联系的优点, 并已经在一些研究中取得了成功。Qu等采用随机森林算法, 重建了南海水下声速数据[15]。Bao 等使用广义回归神经网络 FOAGRNN 模型和果蝇优化算法,估计了太平洋地下盐度数据[16]。Huang等采用了一种基于自主水下航行器(AUV)和水平线阵(HLA)双向通信获得飞行时间并结合人工神经网络的方法, 在浅水区域实现了声速反演[17]。然而, 该方法并未充分考虑水声信道具有严重噪声的特点, 导致其鲁棒性较差不适用于所有区域。此外, 在网络训练模块中, 该方法也没有考虑HLA节点位置间的关联信息对反演性能影响。
针对以上问题和限制, 本文提出了一套针对水声多途信道场景下的声速剖面反演方法, 旨在构建一个具备高精度和鲁棒性的长短期记忆网络(LSTM)模型实现声速剖面反演。为解决传统飞行时间反演方法中的设备同步问题, 本文使用到达时间差替代飞行时间, 并融合到达时间差、信号到达角度(AOA)和节点位置等多模态声学观测量构建感知矩阵, 作为网络输入, 确保模型有效收敛。此外, 本文提出了基于互相关函数(CCF)的硬阈值追踪(HTP)方法, 即CCF-HTP算法, 用于抑制水声信道中的多途效应对时延估计的影响。在此基础上, 本文构建了基于LSTM网络的反演模型, 进一步提升了声速剖面反演的精度和抗噪性能。
1. 场景模型与通信模型
1.1 场景模型
本文的场景模型如图1所示, 主要研究浅海的声速剖面反演问题。在该场景中, 声源被放置在海底, 接收浮标组等距部署在海面附近, 声源发送周期Chirp信号给接收浮标组。接收浮标组主要是由N个配置有矢量水听器和GPS的浮标组成, 所有浮标通过GPS保证同步, 并且构成水平线阵。水平线阵最终通过卫星通信将采集到的数据发送到后端的控制中心进行反演处理。
本文使用由时延信息、角度信息、位置信息构成的感知矩阵实现声速剖面反演, 因此本文的场景需要对这三个要素的获取过程进行定义。
时延信息定义为矩阵
\boldsymbol{ T } , 其由声源发出的Chirp信号到达HLA中每个接收阵元的到达时间差组成的, 该值可通过时延估计方法获得:{\boldsymbol{T}} = \left[ {{t_1} - {t_N},{t_2} - {t_1},{t_3} - {t_2},{t_4} - {t_3}, \cdots ,{t_N} - {t_{N - 1}}} \right], (1) 其中,
N 为HLA阵元数量;{t_n}, \; n = 1,2, \cdots ,N 表示信号到达第n 个HLA阵元之间的信号传播时间。角度信息定义为矩阵
{\boldsymbol{A}} , 其由Chirp信号到达各阵元时的到达角度组成, 该值可通过浮标阵元的矢量水听器测得:{\boldsymbol{A}} = \left[ {{\theta _1},{\theta _2},{\theta _3}, \cdots ,{\theta _N}} \right], (2) 其中,
{\theta _n}, \; n = 1,2, \cdots ,N 表示信号到达第n 个HLA 阵元的角度值。位置信息定义为矩阵
{\boldsymbol{P}} , 其由浮标阵元坐标组成的, 该值可通过GPS获取。为了能够量化场景内各元素的坐标, 需要先建立合适的坐标系。场景对应的坐标系如图1所示, 坐标原点的水平坐标与声源保持一致, 垂直坐标定义为水面并与HLA保持一致, 水平坐标轴从左向右递增, 垂直坐标轴从上向下递增。则{\boldsymbol{P}} 可定义为{\boldsymbol{P}} = \left[ {{M_1},{M_2},{M_3}, \cdots ,{M_N}} \right], (3) 其中,
{M_n}, \; n = 1,2, \cdots ,N 表示第n 个HLA阵元的位置信息, 以其在坐标系里的坐标来表示。1.2 通信模型
由于水声信道存在的多途特性会对到达时间差的获取产生影响, 本文考虑通过建立信道模型的方式来讨论声速剖面反演方法性能。本文采用的信道模型为多途信道模型, 假设发送信号
u(t) 由信源发送至浮标i , 则浮标i 的接收信号{r}_{\text{i}}(t) 可以表示为{r}_{\text{i}}(t) = \sum\limits_{k = 0}^{{D} - 1} {{\alpha _{ik}}} u\left( {t - {\tau _{ik}}} \right) + n(t), (4) 其中,
{D} 为多途数量,{\tau _{ik}} 为声源到浮标i 所经历信道的第k 条路径下的时延,{\alpha _{ik}} 为其对应的幅值衰减,n(t) 为高斯白噪声, 其满足n(t)\sim {{\mathrm{N}}} \left( {0,\sigma _n^2} \right) 。2. 基于长短期记忆网络与射线声学的声速剖面反演方法
2.1 算法整体框架介绍
传统基于飞行时间的反演方法难以实现能耗与精度的平衡, 这是由于水声信道的复杂性使得精确获取飞行时间较为困难。一方面, 由于水声信道存在高延迟大衰减的特性, 通信双方存在同步问题, 通过单向通信获取的飞行时间不可避免地会引入同步误差。另一方面, 较高精度的传播时延往往只能通过双向通信的方式获取, 但这种方式会导致较大的能量损耗, 且仍无法彻底消除同步误差。为了兼顾能耗与精度, 本文提出基于LSTM和射线声学模型的声速剖面反演方法。与传统模型相比, 该方法利用接收节点的到达时间差进行反演, 从而规避传播时延引起的声源与浮标之间的同步问题, 并且仅需单次通信就能完成。此外, 考虑到单模态数据可能无法使模型收敛, 该方法引入接收信号的到达角度、浮标坐标等信息组成感知矩阵进行多模态估计。
该方案的核心思想是利用LSTM网络去挖掘从感知矩阵到声速剖面的非线性映射关系, 尝试建立一个从感知矩阵到声速剖面的通用反演模型。整个反演系统框架如图2所示。系统框架分为三个阶段: 数据生成阶段、训练阶段、测试阶段。数据生成阶段的目的是获取尽可能真实的声速剖面数据, 为此本研究随机选取全球Argo观测系统中记录的历史温度、盐度和压力等数据, 并通过声速经验公式来生成大量的声速剖面训练数据。训练阶段则是对LSTM网络进行训练, 要点在于感知矩阵的获取以及LSTM训练参数的确定。测试阶段则是要验证基于LSTM的反演方法的准确性和可靠性, 在这个阶段本研究将通过数值仿真模拟水下通信过程, 在不同噪声强度的信道影响下测试基于LSTM的反演方案的性能, 为了获取高精度的感知矩阵, 本文还提出了CCF-HTP算法来实现高精度时延估计。由此声速剖面反演工作流程如下:
(1) 使用Argo浮标记录的历史温度、盐度和压力数据, 并通过声速经验公式生成真实的声速剖面数据;
(2) 利用声射线理论计算理论传播时间, 并生成理想的感知矩阵;
(3) LSTM网络训练;
(4) 利用CCF-HTP方法获取测量的感知矩阵;
(5) 利用训练后的LSTM网络进行声速剖面反演进而测试算法性能。
接下来将详细介绍图2中的各个子模块。
2.2 仿真声速剖面数据集生成
全球Argo观测系统由大量布放在全球海洋中自由漂移的Argo剖面浮标组成, 用于测量海洋上层2000 m的温度、盐度和压力等数据。它能对海洋上层的温度、盐度和压力进行连续监测, 其数据可以精确到几小时内。Argo数据虽然不是全海深数据, 但其具有实时性, 可以反映某海域的日变化、季节变化等。利用水温、盐度和压力的Argo观测剖面数据, 结合声速计算的经验公式可以得到声速剖面。部分声速剖面结果如图3所示。
本文主要目的是提出一种泛化能力强的浅海声速剖面反演方法, 并不局限于反演某种类型的声速剖面, 而是适用于绝大多数浅海声速剖面反演的方法。因此本研究选用了从2013年至2022年十年间的全球Argo观测系统(https://argo.ucsd.edu/)中的海域温度、盐度和压力资料数据集, 包含了北太平洋海域各个季节的数据, 并对其中部分典型区域的数据进行分析。
目前国内外常用的声速经验公式有Chen-Millero, Wilson, Leroy以及Frye and Pugh公式, 各个公式所适用的温度、盐度、压力范围有所差异。文献[18]指出, Chen-Millero经验公式[19]在浅海0 m至200 m深度范围内具有较高精度, 因此本文采用该公式计算声速剖面, 具体如下:
\begin{split} c = & {}{\text{ }}{c_W}(T,P) + {c_C}(T,P) + A(T,P)S_\text{al} + \\ & {} B(T,P){S_\text{al}^{3/2}} + D(P){S_\text{al}^2}, \end{split} (5) 式中
\begin{aligned} A(T,P) = &1.389 - 1.262 \times {10^{ - 5}}T + 7.164 \times {10^{ - 5}}{T^2} + 2.006 \times {10^{ - 6}}{T^3} - 3.21 \times {10^{ - 8}}{T^4} +\\ & \left( {9.4742 \times {{10}^{ - 5}}} \right. - 1.2580 \times {10^{ - 5}}T - 6.4885 \times {10^{ - 8}}{T^2} + 1.0507\left. { \times {{10}^{ - 5}}{T^3} - 2.0122 \times {{10}^{ - 10}}{T^4}} \right)P + \\ & \left( { - 3.9064 \times {{10}^{ - 7}}} \right. + 9.1041 \times {10^{ - 9}}T - 1.6002 \times {10^{ - 10}}{T^2} + 7.988\left. { \times {{10}^{ - 12}}{T^3}} \right){P^2} +\\ & \left( {1.100 \times {{10}^{ - 10}} + 6.649 \times {{10}^{ - 12}}T} \right. \left. { - 3.389 \times {{10}^{ - 13}}{T^2}} \right){P^3}, \\ \end{aligned} (6) \begin{aligned} {c_W}(T,P) = &1402.388 + 5.0371T - 5.80852 \times {10^{ - 2}}{T^2} + 3.342 \times {10^{ - 4}}{T^3} - 1.478 \times {10^{ - 6}}{T^4} + \\ & 3.1464 \times {10^{ - 9}}{T^5} + \left( {0.153563 + 6.8982 \times {{10}^{ - 4}}T} \right. - 8.1788 \times {10^{ - 6}}{T^2} + \\ & 1.3621 \times {10^{ - 7}}{T^3}\left. { - 6.1185 \times {{10}^{ - 10}}{T^4}} \right)P + \left( {3.126 \times {{10}^{ - 5}}} \right. - 1.7107 \times {10^{ - 6}}T +\\ & 2.5927 \times {10^{ - 8}}{T^2}\left. { - 2.5335 \times {{10}^{ - 10}}{T^3} + 1.0405 \times {{10}^{ - 12}}{T^4}} \right){P^2}+ \\ & \left( { - 9.7729 \times {{10}^{ - 9}} + 3.8504 \times {{10}^{ - 10}}T} \right.\left. { - 2.3643 \times {{10}^{ - 12}}{T^2}} \right){P^3}, \\ \end{aligned} (7) \begin{split} {c_C}(T,P) =& \left( {0.0029 - 2.19 \times {{10}^{ - 4}}T + 1.4 \times {{10}^{ - 5}}{T^2}} \right)P +\\ & \left( { - 4.76 \times {{10}^{ - 6}} + 3.47 \times {{10}^{ - 7}}T} \right. -\\ & \left. { 2.59 \times {{10}^{ - 8}}{T^2}} \right){P^2} + 2.68 \times {10^{ - 9}}{P^3}, \end{split} (8) \begin{split} B(T,P) = & - 1.922 \times {10^{ - 2}} - 4.42 \times {10^{ - 5}}T + \\ & \left( {7.3637 \times {{10}^{ - 5}} + 1.7945 \times {{10}^{ - 7}}} \right)P, \end{split} (9) D(P) = 1.727 \times {10^{ - 3}} - 7.9836 \times {10^{ - 6}}P, (10) 其中,
T 为温度,^\circ {\text{C}} ;S_\text{al} 为盐度;P 为压力, bar。由于Argo数据并非等间距采样, 因此利用Akima差值将声速剖面数据内插到垂直标准层[20]。部分月份的声速剖面如图3所示, 可以看出本文选择的声速剖面数据包含了各种曲线类型的声速剖面, 使得声速剖面数据集具有泛化性, 从而为后续的声速剖面反演方法提供有效的数据基础。
2.3 基于CCF-HTP的时延估计算法
本文的发送信号为Chirp信号, 其存在明显的自相关特性, 而传统的匹配滤波方法由于多途效应产生的多个相关峰会对时延估计产生干扰, 为了获得更高精度的到达时间差, 本文提出将CCF-HTP算法用于时延估计, 其本质上是一种基于压缩感知(CS)的算法。设发送信号为
u(t) , 接收信号为{u^\prime }(t) , 则互功率谱密度S(w) 可表示为S(f) = \frac{1}{M}S_u^{\text{*}}(f){S_{u'}}(f), (11) 其中,
{S_u}(f) 和{S_{u'}}(f) 分别表示信号u(t) 和{u{'}}(t) 的离散傅里叶变换,M 为离散序列长度。将互功率谱密度S(f) 和互相关函数\xi (\tau ) 分别用矢量形式{\boldsymbol{S}} 和{\boldsymbol{\xi }} 表示, 则可通过维纳–辛钦定理建立压缩感知方程:\boldsymbol{S} = {\boldsymbol{\varPhi \xi }} + \boldsymbol{W}, (12) 其中
\boldsymbol{S} = {[S(0),S({f_1}), \cdots ,S({f_{M - 1}})]^{\text{T}}}, (13) {\boldsymbol{\varPhi }} = \left[ {{\boldsymbol{a}}\left( {{\tau_0}} \right),{\boldsymbol{a}}\left( {{\tau_1}} \right), \cdots ,{\boldsymbol{a}}\left( {{\tau_{N - 1}}} \right)} \right], (14) {\boldsymbol{a}}\left( {{\tau }_{\text{i}}} \right) = {\left[ {1,{{{\mathrm{e}}} ^{{{\mathrm{j}}} 2{\text{π }}{f_1}{\tau }_{\text{i}}}}, \cdots ,{{{\mathrm{e}}} ^{{{\mathrm{j}}} 2{\text{π }}{f_{M - 1}}{\tau }_{\text{i}}}}} \right]^{\text{T}} }, (15) {\boldsymbol{\varPhi }} 表示观测矩阵,\boldsymbol{W} 表示噪声在频域内的噪声矢量,{\tau }_{\text{i}} 是时延差,{f}_{\text{i}} 表示频率。可以证明矩阵
{\boldsymbol{\varPhi }} 满足CS算法要求的RIP性质, 之后可以采用HTP算法求解\xi (\tau ) 。HTP结合了压缩采样匹配追踪算法以及迭代硬阈值算法的优势, 在估计稀疏系数时, 借鉴了压缩采样匹配追踪算法的回溯理念, 能够降低信号的重建误差, 更好地找到支撑位置, 实现稀疏目标信号的重建。根据此性质, 本研究将其运用到压缩感知的时延差估计中, 从而实现更加准确的时延差重建和恢复, 并减少运算时长。CCF-HTP算法流程如表1所示。表 1 CCF-HTP算法流程01 输入: \; \text{iter, }l,u\text{, err, }{u}^{\prime } \begin{array}{l}\text{// iter: }最大迭代次数l\text{: }稀疏度\\ \text{// }u\text{: }发送信号\text{err: }最小误差\text{ }\\ \text{// }{u}^{\prime }\text{: }接收信号\end{array} 02 {\boldsymbol{S}} = \dfrac{1}{M}{\boldsymbol{S}}_{_u}^{\text{H}} \cdot {{\boldsymbol{S}}_{{u^{'}}}} // 计算互功率谱 03 {\boldsymbol{\varPhi }} = \left[ {\boldsymbol{a}\left( {{\tau_0}} \right), \boldsymbol{a}\left( {{\tau_1}} \right), \cdots , \boldsymbol{a}\left( {{\tau_{N - 1}}} \right)} \right] // 计算观测矩阵 04 {{\boldsymbol{r}}^0} = {\boldsymbol{S}}, \, {{\boldsymbol{\xi }}^0} = 0, \, {\mu ^0} = 1 // 初始化参数 05 {\text{while}}(i < {\text{iter}})\,{\text{and}}\,\left( {\left\| {{{\boldsymbol{r}}^n} - {{\boldsymbol{r}}^{n - 1}}} \right\|_{\text{2}}^{\text{2}} > {\text{err}}} \right) 06 \begin{gathered} {{\varGamma }^n} = {\text{\{ indices of }}k{\text{ largest entries of }} \\\qquad {\text{ }}{{\boldsymbol{x}}^n} + {\mu ^n}{{\boldsymbol{\varPhi }}^T}{{\boldsymbol{r}}^{n - 1}}{\text{\} }} \\ \end{gathered} // 取备选更新信号的l个
最大值索引构建索引
集合07 {\boldsymbol{v}} = \arg \min \{ ||{\boldsymbol{S}} - {\boldsymbol{\varPhi v}}||,{\text{supp}} ({\boldsymbol{v}}) \subseteq {{\varGamma }^n}\} // 选择支撑集包含于索引集合的有效信号 08 {{\boldsymbol{\xi }}^n} = {\boldsymbol{v}} // 更新信号 09 {\mu ^n} = \dfrac{{{\text{||(}}{{\boldsymbol{\varPhi }}^{\text{*}}}{\text{(}}{\boldsymbol{S}} - {\boldsymbol{\varPhi }}{{\boldsymbol{x}}^n}{\text{)}}{{\text{)}}_{{{\varGamma }^n}}}{\text{||}}_{\text{2}}^{\text{2}}}}{{{\text{||}}{\boldsymbol{\varPhi }}{{{\text{(}}{{\boldsymbol{\varPhi }}^{\text{*}}}{\text{(}}{\boldsymbol{S}} - {\boldsymbol{\varPhi }}{{\boldsymbol{x}}^n}{\text{))}}}_{{{\varGamma }^n}}}{\text{||}}_{\text{2}}^{\text{2}}}} // 更新步长 10 {{\mathrm{end}}} 11 输出: \; {\boldsymbol\xi }^{n} 2.4 基于射线声学的声速剖面反演机理
本文使用包含到达时间差、到达角度以及阵元坐标的感知矩阵实现声速剖面反演, 其前提是感知矩阵与声速剖面之间存在一定的反演逻辑, 本节通过射线声学模型解释它们之间的联系。
声线跟踪的理论基础为声速分层假设,即任何复杂的声速剖面结构,都可以近似为由多层具有简单结构的声速层组成。此假设是用每个独立层内的声速折线分布来代替整个声速剖面的连续变化。本文使用的声速分层形式为常声速分层, 即每个小层内声速不变, 声线沿直线传播。假设信号的初始入射角度为
{\theta _0} , 初始声速为{c_0} , 由Snell定律和各层声速{c}_{\text{i}} , 可以得到声线在各层中的入射角度有如下关系:\sin {\theta }_{\text{i}} = \frac{{{c}_{\text{i}}}}{{{c_0}}}\sin {\theta _0}, (16) 其中,
{\theta }_{\text{i}} 表示第i 层的入射角度。对于给定的声速剖面, 可以通过等间距采样的方式将其表示为K + 1 个离散声速点的集合, 声线经过完整的K 层后, 到达第n 个接收阵元的水平传播距离计算方法为X_K^n = \sum\limits_{i = 0}^K {\frac{{{c}_{\text{i}}\sin {\theta _0}\left( {{{\textit z}_{i + 1}} - {{\textit z}}_{\text{i}}} \right)}}{{\sqrt {c_0^2 - c_{\text{i}}^2{{\sin }^2}{\theta _0}} }}} , (17) 其中,
{{\textit z}}_{\text{i}} 表示第i 层的深度。经过完整的K层到达第n个接收阵元的传播时间计算方法为T_K^n = \sum\limits_{i = 0}^K {\frac{{{c_0}\left( {{{\textit z}_{i + 1}} - {{\textit z}}_{\text{i}}} \right)}}{{{c}_{\text{i}}\sqrt {c_0^2 - c_{\text{i}}^2{{\sin }^2}{\theta _0}} }}} , (18) 则信号到达第
i 个阵元与第j 个阵元的到达时间差可以表示为\text{TDOA}_{i,j} = T_K^i-T_K^j. (19) 由此可见
\text{TDOA}_{i,j} 、到达角{\theta _0} 以及用于衡量阵元坐标的水平传播距离X_K^n 与由{c}_{\text{i}} 组成的声速剖面之间存在联系, 通过前三者求解后者是可能的, 但这种反演关系难以用闭式解表示, 因此本文使用深度学习方法研究它们之间的联系。2.5 基于LSTM的声速剖面反演方法
本文的场景模型中, 接收端的每个阵元是等间距分布的, 因此信源发出的Chirp直达信号到达每个接收阵元的传播时间存在一定的时序关系, 由此可进一步推断每个直达信号传播时间的差值也存在一定的时序关系。因此, 需要使用一种能够有效处理具有时序关系的输入数据的网络反演声速剖面。
LSTM是特殊的循环神经网络(RNN), 它在RNN基础上增加了遗忘机制和保存机制,可以起到保存较长序列数据中的重要信息, 忽略不重要信息, 以精确的方式传递记忆的作用, 因此十分适合处理时序数据。它最突出的优势就是能够避免长期依赖性问题, 成功解决了原始循环神经网络的缺陷。因此本文使用LSTM网络实现声速剖面反演。
LSTM 在循环神经网络的隐藏层的神经节点中, 增加了一种被称为记忆单元的结构记忆过去的信息, 通过设置遗忘门、输入门和输出门控制历史信息的使用。LSTM单元结构如图4所示,
\boldsymbol{x}_m 为第m 个LSTM单元的输入数据,\boldsymbol{h}_m 为第m 个LSTM单元的输出数据,\boldsymbol{c}_m 为第m 个 LSTM单元的细胞状态。输入门负责决定保留多少当前 LSTM 单元的输入
\boldsymbol {x}_m 到当前的单元状态\boldsymbol{c}_m 中, 计算公式为\boldsymbol{i}_m=\sigma\left(\boldsymbol{W}_{\text{i}}\left[\boldsymbol{h_{m-1}},\boldsymbol{x}_m\right]+\boldsymbol{b}_{\text{i}}\right). (20) 遗忘门决定上一LSTM单元的状态
\boldsymbol{c_{m - 1}} 有多少保留到当前LSTM单元状态\boldsymbol {c}_m 中, 计算公式为\boldsymbol{f}_m = \sigma \left( {\boldsymbol{W}_{\text{f}}\left[ {\boldsymbol{h_{m - 1}},\boldsymbol{x}_m} \right] + \boldsymbol{b}_{\text{f}}} \right), (21) \boldsymbol{\widetilde c}_m = \tanh \left( {\boldsymbol{W_{\text{c}}}\left[ {\boldsymbol{h_{m - 1}},\boldsymbol{x}_m} \right] + \boldsymbol{b_{\text{c}}}} \right), (22) \boldsymbol {c}_m = \boldsymbol{f}_m \odot \boldsymbol{c_{m - 1}} + \boldsymbol{i}_m \odot \boldsymbol{\widetilde c}_m. (23) 输出门可以根据需要有选择地输出信息, 计算公式如下:
\boldsymbol {o}_m =\sigma (\boldsymbol{W}_{\text{o}}[\boldsymbol{h_{m-1}},\boldsymbol{x}_m]+\boldsymbol{b}_{\text{o}}), (24) \boldsymbol{h}_m = \boldsymbol{o}_m \odot \tanh \boldsymbol{c}_m, (25) 式中,
\boldsymbol{W}_{\text{f}},\boldsymbol{W}_{\text{i}},\boldsymbol{W}_{\text{c}},\boldsymbol{W}_{\text{o}} 与\boldsymbol{b}_{\text{f}},\boldsymbol{b}_{\text{i}},\boldsymbol{b}_{\text{c}},\boldsymbol{b}_{\text{o}} 分别为遗忘门、输入门、单元状态、输出门的权重矩阵与偏置项;\boldsymbol{f}_m 为遗忘门输出值;\boldsymbol {\widetilde c}_m 为第m个LSTM单元输入的单元状态;\odot 为哈达玛积运算符;\sigma 为sigmoid激活函数。本文使用的LSTM网络结构如图4所示, 前端是若干个LSTM神经单元, 后端为一个全连接层, 全连接层取最后一个LSTM神经单元的隐藏层输出与声速剖面维度进行映射。
在训练阶段, 将到达时间差、到达角度、阵元位置(POS)数据打包成感知矩阵作为LSTM神经网络的输入, 经过LSTM神经单元处理后输出隐藏层变量, 最后通过一个全连接层映射为声速剖面数据。LSTM根据真实标签声速剖面数据与输出之间的损失值进行训练, 通过不断的迭代训练, 使得LSTM网络学习感知矩阵与声速剖面之间的非线性映射关系, 进而使得生成的声速剖面不断接近真实的声速剖面, 最终训练得到理想的声速剖面反演网络。
在工作阶段, 将测量得到的感知矩阵输入训练后的LSTM神经网络即可得到声速剖面反演结果。由于训练阶段可以离线完成, 且在工作阶段只需要进行一次前向传播计算, 因此大大减少了该方法的计算时间开销。
3. 仿真结果
3.1 仿真平台介绍与仿真参数设置
本文使用数值仿真验证了算法性能, 所有的实验均在硬件为Intel(R) Xeon(R) Silver 4210 CPU @ 2.20 GHz, 256 GB内存, NVIDIA GeForce RTX 3090, CUDA 11.7; 软件为Python 3.9的环境中进行。仿真实验在深度200 m, 水平距离为200 m的范围内进行。 声源位置为(x = 0 m, y = 200 m), 固定于海底, 接收浮标的第一个位置为(x = 10 m, y = 0 m), 最后一个位置为(x = 200 m, y = 0 m), 相邻浮标间距10 m固定排列, 且假设浮标坐标位置固定不变。多途信号根据前文所述内容由BELLHOP模型生成, 并在不同信噪比下进行实验, 模型参数设置如表2所示。
表 2 参数设置参数 仿真取值 采样率 (kHz) 100 Chirp信号频率范围 (kHz) 20~30 Chirp信号周期 (ms) 40 水声信道多途数量 10 最大时延拓展 (ms) 30 输入层神经元个数 4 输出层神经元个数 25 接收浮标个数 20 PSO粒子数 50 PSO迭代次数 50 学习率 0.0001 声速剖面训练集数量 8800 声速剖面测试集数量 2200 声速剖面深度采样间隔 (m) 8 训练迭代次数 10000 3.2 时延估计算法性能验证
为了验证CCF-HTP算法在时延估计方面的有效性,通过仿真实验评估该算法在不同信噪比影响下的时延估计准确度。仿真发送的信号为Chirp信号, 其具体参数如表2所示, 其时频特性与自相关特性如图5所示。
仿真实验随机取10组不同的声速剖面分布, 信源与每个浮标之间的信道根据BELLHOP信道模型生成, 并在不同信噪比下根据式(4)计算得到接收信号。图6为当信噪比取值−5 dB时的接收信号的时频特性与自相关特性图。
之后利用TDE算法处理接收信号, 即可获取时延估计采样点数, 利用采样点数误差
{E_{{\text{SP}}}} 评价算法时延估计性能:{E_{{\text{SP}}}} = \left| {{N_{{\text{TP}}}} - {N_{{\text{TR}}}}} \right|, (26) 其中,
{N_{{\text{TP}}}} 为传播时延估计值,{N_{{\text{TR}}}} 为传播时延真实值。为了进一步验证算法性能, 将CCF-HTP算法与基于正交匹配追踪(OMP)以及匹配滤波(MF)取阈值的时延估计(TDE)算法在不同信噪比下进行对比, 每个信噪比下都重复100次实验后取均值, 仿真结果如图7(a)所示。较高信噪比时三种TDE方法的表现相差无几, 这是因为高信噪比下的相关峰较容易分辨; 而低信噪比时传统算法的性能急剧下降, 但CCF-HTP算法有较高的估计精度和较好的稳定性。由此可见, 相较于传统的MF和OMP算法, CCF-HTP算法在时延估计方面具备更好的性能。
为了能够更加直观地呈现CCF-HTP算法的抗多途特性, 仿真实验在不同距离以及不同多途数量的环境下进行TDE计算, 其中声波的传播损失定义为球面波的扩展损失与吸收损失之和[21], 并由此确定加性噪声, 假设初始距离为200 m且对应的信噪比为15 dB, 每种距离和多途数量的组合都进行50次实验取均值, 其仿真结果如图7(b)所示。从图中可以看出, 无论是当声波传播距离还是多途数量增加时, 本文所提算法的性能均优于传统算法。
3.3 反演算法性能验证
本文通过仿真实验验证LSTM的反演性能, 相关的实验平台和参数设置如3.1节所述。所使用的数据集共包含了11000条声速剖面曲线信息, 深度范围为200 m, 采样点为25。其中, 随机选取了8800个数据作为训练集, 2200个数据作为测试集。
该部分的仿真实验主要分为训练阶段和测试阶段。训练阶段旨在优化LSTM网络的学习性能, 因此需要设置合适的参数。为了获取最优参数, 该阶段在不同的LSTM单元数和隐藏层神经单元数组合下训练LSTM网络, 并选择使损失值最小的组合作为LSTM网络的最终参数。该实验的训练迭代次数统一设置为10000次, 其他参数如表2所示, 并采用均方根误差(RMSE)函数作为损失函数来衡量损失值, 仿真结果如图8所示。
图8(a)是不同LSTM单元数和隐藏层神经单元数组合下的网络训练损失结果, 当隐藏层个数取100, LSTM单元个数取4时损失值最小。因此, 后续测试将采用该参数组合下的LSTM模型进行验证。图8(b)是在该参数组合下, LSTM网络在训练集和测试集上损失值随训练迭代次数增加的变化情况。从图中可知, 训练集和测试集的损失值都随着迭代次数增加而不断减小, 最终趋于稳定。这表明: 一方面, 测试集和训练集的最小平均损失值小于2 m/s, 可以认为模型达到收敛; 另一方面, 该参数组合下的LSTM网络的训练阶段并没有出现欠拟合或者过拟合的现象, 因此该模型可以用于后续的测试阶段。
为了验证本文所提的使用多模态感知矩阵的重要性, 仿真通过改变模型的输入参数研究不同输入参数对训练过程的影响。输入参数分别调整为(TDOA, POS)、(AOA, POS)以及感知矩阵并在相同的数据集下进行训练, 迭代次数固定为10000次, 结果如图9所示。从图中可以看出, 单独使用到达时间差或者到达角度作为输入参数LSTM模型都无法收敛, 只有将到达角度与到达时间差组合而成的感知矩阵作为输入时, LSTM模型才能学习到输入与输出之间的非线性关系并收敛。这一结果证明了本文所用的多模态感知矩阵的必要性。
测试阶段旨在评估LSTM网络的反演精度和抗噪性能, 该阶段将理想感知矩阵(由原始测试集数据生成)以及测量感知矩阵(由测试集经过3.2节的TDE过程处理后生成)分别输入训练好的LSTM网络中进行反演, 并将其结果与测试集的声速剖面标签进行对比, 衡量反演精度的指标为RMSE。为了进一步证明算法的有效性, 本文使用粒子群优化(PSO)算法用于对比。图10是LSTM网络对理想感知矩阵和测量感知矩阵的反演结果, 图中的两条IDEA曲线表示的是理想感知矩阵的处理结果, 因此与信噪比无关, 可以看到LSTM的反演精度要好过传统的PSO算法, 图中的两条MEAS曲线表示的是测量感知矩阵的处理结果, 可以看出, LSTM在处理带噪数据时的效果优于PSO算法, 这表明LSTM具有较好的抗噪能力。
图11是输入为理想感知矩阵时的部分反演结果对比, 从图中可以看出LSTM模型的反演结果比传统PSO优化算法的反演结果要更加贴近于原始剖面。这两个样本的反演误差如图12(b)(c)所示, 与传统方法相比, 基于LSTM的反演剖面的误差方差要更小。图12(a)为LSTM网络处理理想感知矩阵时测试集数据的误差分布直方图, 从图中可以看出, 大部分测试数据的反演误差小于4 m/s, 因此可认为LSTM网络模型具备一定的泛化能力, 可以实现大部分声速剖面数据的反演。
4. 结论
本文提出了基于LSTM网络和射线声学的区域海洋声速剖面反演方法并通过数值模拟的方式对反演性能进行了探讨。仿真结果表明，在不同声波传播距离与多途数量的时延估计中， CCF-HTP时延估计算法抗多途性能优于传统的时延估计算法；在不同输入参数组合的声速剖面反演中，只有当同时使用到达时间差和到达角度作为输入参数时，LSTM 模型才能够有效收敛，这说明多模态参数对于此模型的准确性和收敛性至关重要，并且与传统优化算法相比，基于LSTM的声速剖面反演方法在水声多途场景中表现出更好的鲁棒性。
-
表 1 CCF-HTP算法流程
01 输入: \; \text{iter, }l,u\text{, err, }{u}^{\prime } \begin{array}{l}\text{// iter: }最大迭代次数l\text{: }稀疏度\\ \text{// }u\text{: }发送信号\text{err: }最小误差\text{ }\\ \text{// }{u}^{\prime }\text{: }接收信号\end{array} 02 {\boldsymbol{S}} = \dfrac{1}{M}{\boldsymbol{S}}_{_u}^{\text{H}} \cdot {{\boldsymbol{S}}_{{u^{'}}}} // 计算互功率谱 03 {\boldsymbol{\varPhi }} = \left[ {\boldsymbol{a}\left( {{\tau_0}} \right), \boldsymbol{a}\left( {{\tau_1}} \right), \cdots , \boldsymbol{a}\left( {{\tau_{N - 1}}} \right)} \right] // 计算观测矩阵 04 {{\boldsymbol{r}}^0} = {\boldsymbol{S}}, \, {{\boldsymbol{\xi }}^0} = 0, \, {\mu ^0} = 1 // 初始化参数 05 {\text{while}}(i < {\text{iter}})\,{\text{and}}\,\left( {\left\| {{{\boldsymbol{r}}^n} - {{\boldsymbol{r}}^{n - 1}}} \right\|_{\text{2}}^{\text{2}} > {\text{err}}} \right) 06 \begin{gathered} {{\varGamma }^n} = {\text{\{ indices of }}k{\text{ largest entries of }} \\\qquad {\text{ }}{{\boldsymbol{x}}^n} + {\mu ^n}{{\boldsymbol{\varPhi }}^T}{{\boldsymbol{r}}^{n - 1}}{\text{\} }} \\ \end{gathered} // 取备选更新信号的l个
最大值索引构建索引
集合07 {\boldsymbol{v}} = \arg \min \{ ||{\boldsymbol{S}} - {\boldsymbol{\varPhi v}}||,{\text{supp}} ({\boldsymbol{v}}) \subseteq {{\varGamma }^n}\} // 选择支撑集包含于索引集合的有效信号 08 {{\boldsymbol{\xi }}^n} = {\boldsymbol{v}} // 更新信号 09 {\mu ^n} = \dfrac{{{\text{||(}}{{\boldsymbol{\varPhi }}^{\text{*}}}{\text{(}}{\boldsymbol{S}} - {\boldsymbol{\varPhi }}{{\boldsymbol{x}}^n}{\text{)}}{{\text{)}}_{{{\varGamma }^n}}}{\text{||}}_{\text{2}}^{\text{2}}}}{{{\text{||}}{\boldsymbol{\varPhi }}{{{\text{(}}{{\boldsymbol{\varPhi }}^{\text{*}}}{\text{(}}{\boldsymbol{S}} - {\boldsymbol{\varPhi }}{{\boldsymbol{x}}^n}{\text{))}}}_{{{\varGamma }^n}}}{\text{||}}_{\text{2}}^{\text{2}}}} // 更新步长 10 {{\mathrm{end}}} 11 输出: \; {\boldsymbol\xi }^{n} 
下载: 导出CSV
表 2 参数设置
参数 仿真取值 采样率 (kHz) 100 Chirp信号频率范围 (kHz) 20~30 Chirp信号周期 (ms) 40 水声信道多途数量 10 最大时延拓展 (ms) 30 输入层神经元个数 4 输出层神经元个数 25 接收浮标个数 20 PSO粒子数 50 PSO迭代次数 50 学习率 0.0001 声速剖面训练集数量 8800 声速剖面测试集数量 2200 声速剖面深度采样间隔 (m) 8 训练迭代次数 10000
下载: 导出CSV
-
[1] Fuda J L, Millot C, Taupier-Letage I, et al. XBT monitoring of a meridian section across the western Mediterranean Sea. Deep Sea Res. Part I, 2000; 47(11): 2191−2218 DOI: 10.1016/S0967-0637(00)00018-2
[2] Gettelman A, Geer A J, Forbes R M, et al. The future of Earth system prediction: advances in model-data fusion. Sci. Adv., 2022; 8(14): eabn3488 DOI: 10.1126/sciadv.abn3488
[3] 张忠兵, 马远良, 杨坤德, 等. 浅海声速剖面的匹配波束反演方法. 声学学报, 2005; 30(2): 103−107 DOI: 10.3321/j.issn:0371-0025.2005.02.002 [4] Tolstoy A, Sotirin B. Acoustic tomography via matched field processing. J. Acoust. Soc. Am, 1995; 89(5): 393−406 DOI: 10.1121/1.411711
[5] Shang E C. Ocean acoustic tomography based on adiabatic mode theory. J. Acoust. Soc. Am, 1989; 85(4): 1531−1537 DOI: 10.1121/1.397355
[6] Roux P, Iturbe I, Nicolas B, et al. Travel-time tomography in shallow water: Experimental demonstration at an ultrasonic scale. J. Acoust. Soc. Am, 2011; 130(3): 1232−1241 DOI: 10.1121/1.3621271
[7] Li F H, Zhang R H. Inversion for sound speed profile by using a bottom mounted horizontal line array in shallow eater. Chin. Phys. Lett., 2010; 27(8): 084303 DOI: 10.1088/0256-307X/27/8/084303
[8] Voronovich A G, Shang E C. Numerical simulations with horizontal-refraction-modal tomography. J. Acoust. Soc. Am, 1997; 101(5): 2636−2643 DOI: 10.1121/1.418504
[9] Skarsoulis E K, Athanassoulis G A, Send U. Ocean acoustic tomography based on peak arrivals. J. Acoust. Soc. Am, 1996; 100(2): 797−813 DOI: 10.1121/1.416212
[10] Yang T C, Yates T. Matched-beam processing: Application to a horizontal line array in shallow water. J. Acoust. Soc. Am, 1998; 104(3): 1316−1330 DOI: 10.1121/1.424341
[11] 廖光洪, 朱小华, 林巨, 等. 海洋声层析观测技术和方法. 海洋学报, 2010; 32(3): 14−22 [12] 唐俊峰, 杨士莪. 由传播时间反演海水中的声速剖面. 哈尔滨工程大学学报, 2006; 27(5): 733−736 DOI: 10.3969/j.issn.1006-7043.2006.05.022 [13] 李鹏程, 冯海泓, 李记龙. 基于参考点之间水声传播的声速剖面反演. 声学技术, 2023; 42(4): 446−451 DOI: 10.16300/j.cnki.1000-3630.2023.04.006 [14] 张忠兵, 马远良, 倪晋平. 基于声线到达时差的浅海声速剖面反演. 西北工业大学学报, 2002; 20(1): 36−39 DOI: 10.3969/j.issn.1000-2758.2002.01.009 [15] Liu C, Qu K. Wide-area sound speed profile estimation based on a pre-classification scheme for sound speed perturbation modes. Front. Marine Sci., 2023; 10: 1130061 DOI: 10.3389/fmars.2023.1130061
[16] Chen Z, Wang P, Bao S, et al. Rapid reconstruction of temperature and salinity fields based on machine learning and the assimilation application. Front Marine Sci., 2022; 9: 985048 DOI: 10.3389/fmars.2022.985048
[17] Huang W, Li D, Jiang P. Underwater sound speed inversion by joint artificial neural network and ray theory. Proceedings of the 13th International Conference on Underwater Networks & Systems, ACM, Shenzhen, China, 2018: 1−8
[18] 吴碧, 陈长安, 林龙. 多种声速经验公式的比较分析研究. 全国水声学学术会议, 中国声学学会, 广东湛江, 2013: 300−302 [19] 王子蘅, 王振杰, 聂志喜, 等. 声速剖面EOF重构的实测数据采样深度研究. 海洋科学, 2021; 45(6): 126−134 DOI: 10.11759/hykx20200711001 [20] Xia H, Yang K, Duan R, et al. Analysis and estimation of critical depth based on SODA datasets in the Philippine Sea. IEEE/OES China Ocean Acoustics, IEEE, Harbin, China, 2016: 1−7
[21] 刘伯胜. 水声学原理. 北京: 科学出版社, 2019
计量
- 文章访问数: 220
- HTML全文浏览量: 31
- PDF下载量: 98
