海水中的声速决定了海洋中的声传播特性^[1]，剖面受到很多海水环境参数的影响，其中有一些相对稳定的参数，如海底的结构、海洋的深度等；还有一些动态变化的参数，如海水的生物群、海流、温度、盐度等。这些动态变化的参数，会随着时间的变化而发生变化。变化不稳定，会导致海水中声速产生变化，使得水层声速剖面具有明显的时空变化。

海水声速剖面的显著时间演化特性为声速剖面的序贯反演提供了可行性。Candy^[2]将扩展卡尔曼滤波算法引入到海洋声学反演中，将声速剖面的时间演化过程建立为一阶时间演化模型。Yardim^[3]通过经验正交函数(empirical orthogonal function，EOF)描述声速剖面，采用序贯滤波中的卡尔曼滤波和粒子滤波的方法对地声参数进行跟踪。Carriere^[4-5]在水平变化环境中将声学反演问题改进为状态-空间模型，通过序贯滤波的数据同化算法进行海洋环境状态估计。李建龙等^[6]通过海洋动力学模型与声速的关系反演声速场的时间与空间结构信息并验证了算法的有效性；金丽玲等^[7]利用自回归分析拟合方法将声速场扰动建模为高阶自回归演化模型，通过集合卡尔曼滤波序贯的估计时间演化的海洋声速场。这些研究将时变声速剖面的反演问题改进为由描述声速剖面时变特性的状态方程与包含测量声场信息的测量方程组成的状态-空间模型，应对声场环境的时变特性。

由于浅海声场环境的复杂性，难以直接描述声速剖面时变特性的状态方程，通常为一阶随机游走过程，此时的状态方程不能很好的预测状态变量的走向，序贯滤波的应用会受到限制。而深度学习的分析归纳能力强，可以将复杂的问题表示为嵌套的层次概念体系，所以尝试利用深度学习方法的循环神经网络，处理相关数据，描述声速剖面的时变特性。学习数据构建模型并配合序贯滤波器进行估计，这一方法在其他学科中广泛应用。例如，文献[8]中将长短期记忆人工神经网络(long short term memory network，LSTM)与卡尔曼滤波(Kalman filter，KF)深度融合，提出了一种KF修正的LSTM混合模型KF-LSTM，其结果与LSTM模型、浅层神经网络模型和KF模型单独使用的结果相比，能够有效提升性能；文献[9]中将BP神经网络用于测量方程的优化，再通过卡尔曼滤波器进行估计；文献[10]利用神经网络对非线性系统建立状态空间模型，采用高阶容积卡尔曼滤波对新的状态进行实时更新，从而达到神经网络对非线性系统模型的真实逼近以及对状态值的精确估计。

本文将深度学习中的循环神经网络引入到序贯滤波中，利用循环神经网络学习历史水文数据，对浅海环境下的时变声速剖面进行学习建模，在此基础上，利用集合卡尔曼滤波进行了对声速剖面和声源位置的联合反演，通过实测声速剖面的仿真声学数据验证了方法的可行性。

1 声速剖面时变特性模型

LSTM神经网络由Hochreiter等^[11]提出，并由Graves^[12]进行改进，是基于RNN的一种完善，解决RNN中易出现的梯度消亡问题。在LSTM中，通过门结构来对细胞状态信息进行增加或删除。门结构通常由一个sigmoid神经网络层和逐点乘积的操作组成，从而选择性地让信息通过。LSTM通过输入门、遗忘门和输出门3个门结构来实现对信息的控制。

基于循环神经网络进行模型构建与训练。模型训练的数据集来自于某次浅海实验中A、B 2点约10 d的水文数据，数据集中包括A点声速剖面和B点声速剖面，每10 s记录一次数据。利用这个数据集构建浅海环境下的声速剖面时变特性模型，最终达到已知A点处前1 h的水文数据，可以预测B点1 h后的声速剖面。在Keras框架下进行了模型的构建与训练。

在模型训练结束后，使用测试集数据进行模型测试。图 1为实际声速剖面、模型预测的声速剖面、模型预测的误差。可以看出，模型预测值的变化趋势，与实际声速剖面的时间演化趋势基本一致，只在温跃层深度处，因声速剖面起伏较大，模型的预测误差变大。图 2为测试集数据进行预测的均方根误差，在温跃层起伏较大时，均方根误差较大。

2 基于声速剖面时变特性模型的联合反演 2.1 声速剖面参数化

EOF方法能够有效地将声速剖面参数化，降低表示声速剖面所需要参数个数^[13-14]。假设在一段时间内测量得到N条声速剖面SSP数据，在深度上进行等间隔线性插值。离散到M个深度上：c₁(z_j), c₂(z_j), …, c_N(z_j), j=1, 2, …, M

从而可以通过求平均值的办法得到平均声速剖面：

$ \bar{c}(z)=\frac{1}{N} \sum\limits_{i=1}^{N} c_{i}(z) $

(1)

则第i条SSP在第j个深度处的声速起伏为：

$ \Delta c_{i}\left(z_{j}\right)=c_{i}\left(z_{j}\right)-\bar{c}\left(z_{j}\right), i=1,2, \cdots, N $

(2)

协方差矩阵R可以表示为：

$ \boldsymbol{R}=\frac{1}{N} \boldsymbol{C} \boldsymbol{C}^{\mathrm{T}} $

(3)

$ \boldsymbol{C}=\left[\begin{array}{cccc} \Delta c_{1}\left(z_{1}\right) & \Delta c_{2}\left(z_{1}\right) & \cdots & \Delta c_{N}\left(z_{1}\right) \\ \Delta c_{1}\left(z_{2}\right) & \Delta c_{2}\left(z_{2}\right) & \cdots & \Delta c_{N}\left(z_{2}\right) \\ \vdots & \vdots & \vdots & \vdots \\ \Delta c_{1}\left(z_{M}\right) & \Delta c_{2}\left(z_{M}\right) & \cdots & \Delta c_{N}\left(z_{M}\right) \end{array}\right] $

(4)

对矩阵R进行特征值分解，可以得到：

$ \boldsymbol{R}=\sum\limits_{n=1}^{M} \lambda_{n} \boldsymbol{f}_{n} \boldsymbol{f}_{n}^{\mathrm{T}} $

(5)

式中λ₁>λ₂>…>λ_M为对应特征向量f_n的特征值。

由于本征值快速衰减，通常前几阶数EOF系数就可以较为准确地描述声速剖面：

$ \boldsymbol{c}(z)=\overline{\boldsymbol{c}}(z)+\sum\limits_{k=1}^{K} \boldsymbol{\alpha}_{k} \boldsymbol{f}_{k}(z) $

(6)

式中α_k和f_k(z)分别为第k阶EOF系数和第k阶EOF向量。

2.2 空间-状态方程

本文利用集合卡尔曼滤波方法^[15]对浅海环境下时变声速剖面及声源位置同时进行反演，系统的状态方程需要同时描述声速剖面和声源位置的变化。在先验数据的基础上，利用前3阶经验正交系数表示声速剖面:

$ \boldsymbol{\alpha}_{k \mid k-1} \equiv\left[\begin{array}{l} \alpha_{1} \\ \alpha_{2} \\ \alpha_{3} \end{array}\right]_{k} \equiv L\left(\boldsymbol{\alpha}_{k-1}\right)+\boldsymbol{w}_{k} $

(7)

式中: w_k为状态噪声向量；L表示之前所构建的声速剖面时变特性模型，以A点前1 h的声速剖面数据作为先验数据，逐步预测B点声速剖面，并进行EOF分解，更新α。

$ \boldsymbol{s}_{k \mid k-1} \equiv\left[\begin{array}{c} Z_{k \mid k-1} \\ R_{k \mid k-1} \\ v_{k 1 k-1} \end{array}\right] \equiv \boldsymbol{F}\left(\boldsymbol{s}_{k-1}\right)+\boldsymbol{V}\left[\begin{array}{c} v_{d} \\ v_{a} \end{array}\right], k=1,2, \cdots, K $

(8)

$ \boldsymbol{F}=\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & \Delta t \\ 0 & 0 & 1 \end{array}\right], \boldsymbol{V}=\left[\begin{array}{cc} 1 & 0 \\ 0 & \Delta t^{2} / 2 \\ 0 & \Delta t \end{array}\right] $

(9)

式中：s为声源的参数，Z为声源的深度，R为水平距离；v为速度；Δt为2次测量之间的时间差；v_d和v_a分别是声源深度与速度的随机变量。

则系统的状态向量定义及系统的状态方程和测量方程为：

$ \boldsymbol{x}_{k}^{\mathrm{T}}=\left[\begin{array}{cc} \boldsymbol{s}_{k}^{\mathrm{T}} & \boldsymbol{\alpha}_{k}^{\mathrm{T}} \end{array}\right] $

(10)

$ \boldsymbol{x}_{k}=f\left(\boldsymbol{x}_{k-1}\right)+\boldsymbol{v}_{k} $

(11)

$ \boldsymbol{y}_{k}=h\left(\boldsymbol{x}_{k}\right)+\boldsymbol{w}_{k} $

(12)

式中：f(·)为状态向量的状态转移函数；L为声速剖面时变特性模型；F为声源位置状态转移矩阵。v_k和w_k为过程噪声和测量噪声；h(·)为简正波模型；考虑远场近似，k时刻水平距离R、深度Z处的声源，在接受阵列处的声压P可以表示为：

$ P(R, z)=\frac{\mathrm{i}}{\sqrt{8 {\rm{ \mathsf{ π} }} \rho}} \mathrm{e}^{\frac{\mathrm{i} {\rm{ \mathit{ π} }}}{4}} \sum\limits_{m=1}^{M} \varPhi_{m}(Z) \varPhi_{m}(z) \frac{\mathrm{e}^{\mathrm{i} k_{m} R}}{\sqrt{k_{r m}}} $

(13)

式中：Φ_m(z)为本征函数，描述了第m号简正波的模式形状；k_rm为本征值，即水平波数；ρ为海水密度；M为简正波个数；Z为声源深度。

2.3 集合卡尔曼滤波

集合卡尔曼滤波^[16]为基于蒙特卡罗采样的卡尔曼滤波，可以解决非线性和非高斯的问题。本文在2.2节的基础上，获得由N个样本点组成的背景集合X_k和观测集合Z_k，并计算背景集合的平均状态分别为：

$ \boldsymbol{X}_{k}=\left\{\boldsymbol{x}_{k, i}, i=1,2, \cdots, N\right\} $

(14)

$ \boldsymbol{Z}_{k}=\left\{\boldsymbol{z}_{k, i}, i=1,2, \cdots, N\right\} $

(15)

$ \boldsymbol{x}_{k, i}=f\left(\boldsymbol{x}_{k-1}\right)+\boldsymbol{v}_{k} $

(16)

$ \bar{\boldsymbol{x}}_{k}=\frac{1}{N} \sum\limits_{i=1}^{N} \boldsymbol{x}_{k, i} $

(17)

其中观测集合在背景集合的基础上构建：

$ \boldsymbol{z}_{k, i}=h\left(\boldsymbol{x}_{k, i}\right) $

(18)

并通过平均的方式获得状态误差的协方差矩阵P_xh^k和测量误差的协方差矩阵P_hh^k。

利用最新的观测信息，通过卡尔曼滤波更新当前时刻的背景集合里的每个样本点，得到分析集合：

$ \boldsymbol{X}_{k}^{a}=\left\{\boldsymbol{x}_{k, i}^{a}, i=1,2, \cdots, N\right\} $

(19)

$ \boldsymbol{x}_{k, i}^{a}=\boldsymbol{x}_{k, i}+\boldsymbol{G}_{k}\left(\boldsymbol{z}_{k, i}-h\left(\boldsymbol{x}_{k, i}\right)\right) $

(20)

其中G_k为卡尔曼增益：

$ \boldsymbol{G}_{k}=\boldsymbol{P}_{x h}^{k}\left(\boldsymbol{P}_{h h}^{k}\right)^{-1} $

(21)

利用分析集合的平均值对k时刻的状态进行更新：

$ \boldsymbol{x}_{k}=\frac{1}{N} \sum\limits_{i=1}^{N} \boldsymbol{x}_{k, i}^{a} $

(22)

3 联合反演结果 3.1 匀速直线运动声源仿真

假设声源做匀速直线运动，利用KRAKEN声场计算工具，使用实测声速剖面数据进行浅海波导下的声场仿真。仿真波导的海深为84 m，海底密度为1.78 g/cm³，海底声速为1 650 m/s，海底衰减为0.13 dB/λ。

利用LSTM声速剖面时变特性模型进行了声速剖面和声源位置的联合反演，记为LSTM-EnKF方法；同时利用了一阶时间演化模型进行了声速剖面和声源位置的联合反演，记为EnKF方法。将反演结果进行对比。

LSTM-EnKF方法和EnKF方法下声源参数的反演结果如图 3、图 4所示。图 3(a)、(c)、(e)分别为LSTM-EnKF方法下声源深度、距离、速度的反演结果与实际的对比，图 3(b)、(d)、(f)为3个参数的反演误差；图 4(a)、(c)、(e)分别为EnKF方法下声源深度、距离、速度的反演结果与实际的对比，图 4(b)、(d)、(f)为3个参数的反演误差。可以发现2种方法都能够对声源参数进行追踪，且均展现了较好的追踪性能，3个参数的误差均在0左右。这是由于在进行声源追踪时，2种方法使用了相同的状态方程和状态转移矩阵。

如图 5所示，对于EOF前3阶系数的追踪结果，在第3阶EOF系数的追踪上，EnKF追踪结果的误差要略大于LSTM-EnKF追踪结果的误差。并且可以看出，使用一阶时间演化模型的EnKF方法，在前30 min左右出现了较大的误差，也正是EnKF方法在第3阶EOF系数追踪上，误差较大的时刻。此外，图 4显示EOF3的跟踪比EOF1和EOF2系数差。与前2阶相比，第3阶EOF系数对SSP扰动的控制作用较小，导致声速的微小变化。通过比较，验证了LSTM-EnKF法在声速剖面的追踪上的可行性，LSTM-EnKF法能够估计所有参数，且波动范围较小，与真实值几乎完全吻合。

3.2 非匀速运动声源仿真

在相同仿真参数下，利用KRAKEN声场计算工具，使用实测声速剖面数据，选取该海域某次浅海实验移动声源的实测运动参数，进行非匀速运动声源的浅海声场仿真。使用LSTM-EnKF方法，进行了声速剖面和声源位置的联合反演。如图 7所示，与图 3中LSTM-EnKF方法对匀速直线运动声源的追踪结果相比较，对于实测声源的追踪，由于声源做非匀速运动，追踪结果的误差略大。其中声源速度的追踪性能较差，这可能是由于在表1所示的仿真参数下，测量方程(18)所对应的代价函数对声源速度的敏感度较低。

LSTM-EnKF方法下声源参数的反演结果如图 8所示, 上中下3个参数分别为声源深度、水平距离及声源速度，其中进行了估计值与真实值的比较，并展示了声源参数的反演误差。图 8、9中为EOF 3阶系数与声速剖面的反演结果，第3阶EOF系数追踪效果较差，但由于第3阶EOF系数对SSP扰动的控制作用较小，声速剖面的时间演化过程仍能够很好地描述；其误差与图 6(b)相比略小，对比图 5(b)与图 8，这可能是由于第2阶EOF系数对SSP扰动的控制要大于第3阶系数。

4 结论

1) 使用循环神经网络对历史水文数据进行学习建模，能够得到对声速剖面的时变特性较为精确的描述。

2) 在相同的仿真数据下，使用循环神经网络训练得到的模型作为状态方程，相较于使用一阶随机游走过程，前者联合反演结果更为精确，能够追踪声速剖面波动。

在仿真条件中，声源的移动轨迹较为平缓，应使用复杂多变的运动轨迹来进行方法的验证；仿真过程和联合反演过程中将浅海环境视为水平均匀环境，未考虑声速剖面空间变化对联合反演结果的影响；该方法需要测量点前一小时的声速剖面数据作为先验数据。在后续的工作中，要根据以上几点，探寻利用深度学习，如何进行水平非均匀环境下时变的声速剖面的描述，并进行相应复杂环境下的目标定位。