基于内波统计特性的声速剖面预测方法

2021-06-10苏林孙炳文胡涛任群言王文博郭圣明马力

哈尔滨工程大学学报 2021年6期

苏林，孙炳文，胡涛，任群言，王文博，郭圣明，马力

(1.中国科学院声学研究所中国科学院水声环境特性重点实验室，北京 100190； 2.中国科学院大学，北京 100049)

海洋中的声速剖面是影响水声装备性能的重要因素之一，决定了海洋中的声传播特性。海洋内波是广泛存在于近海大陆架的一种中尺度海洋现象，在夏季，由于海水沿重力方向的温度分布不均匀，海水中会形成一定的密度分层[1]。内波活动使海水的声速分布产生剧烈变化，严重影响水下声场分布[2-3]，进而干扰水下声场预报效能[4]，特别是夏季有温跃层时，若声速扰动改变了跃层的斜率，可使声波向下折射到海底进而影响探测作用距离[5]，因此声速剖面的实时监测非常重要，为了预先了解区域声场分布，对声速剖面的预报更有实用价值。

海水声速与温度、盐度和静压力有关，随温度、盐度、压力的增加而增加，其中以温度的影响最为显著，同时，温度的变化也是最为显著的，通常可通过现场监测的温度数据和声速经验公式计算得到声速剖面分布。声速剖面具有明显的深度分布，常表示为深度和声速的矩阵数组，涉及参数量较大，即便在浅海环境也给声速剖面的预报带来很多困难，为简化表达，提出了经验正交函数(EOF)的方法[6]，基函数确定后，用有限的几阶系数即可重构声速剖面，可以大大降低所需要的参数。Gaelle[7]在对意大利墨西拿海峡区域的内波分布研究中，指出可以利用2个EOF系数散点分布图的统计特性监测内波，利用该特征可以进一步的降低声速剖面预报所需参数，提升预报效率。

目前海温预报主要方式有经验预报方法、统计预报方法[8]和数值预报方法。其中经验预报方法与人的主观因素有很大关系；数值预报运算量极大，在要求预报时效性的基础上对计算环境要求较高，且在初始场构建及物理过程参数化等方面有待完善；统计方法则是利用以往历史数据或在线的实时、准实时数据，通过概率统计的方法计算预报未来海水温度的变化规律。卡尔曼滤波方法是以统计理论为基础的代表性同化方法，现在已逐渐应用到海洋数值模拟和预报中[9]，但卡尔曼滤波方法在状态空间建模时仍需要较为准确的状态方程描述海洋变化过程。神经网络具有较好的自学习能力，可跳过复杂的海洋模型，直接从历史及在线数据中挖掘内在联系，尤其近年来海洋预报数据量骤增，积累了大量的潜浮标监测数据，均为神经网络在海洋预报中的应用[10-11]提供了有力支撑。长短时记忆循环神经网络(LSTM)[12]能有效利用时间序列数据中长距离依赖信息的能力，非常适合海洋温度、声速等时间序列的预测[13]。

本文深入挖掘了有无内波情形下时序声速剖面所呈现的不同统计特征，结合LSTM对长时间序列的记忆功能，通过深度学习中神经网络的方法学习历史水文数据的变化规律，对浅海典型水文现象影响下的声速剖面变化进行实时预测。

1 浅海内波环境下的声速剖面统计特征分析

1.1 南海试验概况

2019年9月，中科院声学所水声环境特性重点实验室在中国南海海域组织了海洋调查试验，试验海区地形如图1所示，海深约85 m，试验期间在沿等深线方向布放了2套温度链潜标，西侧站位记为B站位，东侧站位记为O站位。

图1 南海试验地形Fig.1 Bathymetry of the experiment area in South China Sea

试验从9月8日持续到9月15日，共采集到约171 h的温度数据，如图2所示，在整个试验周期内可以观测到丰富的内波活动，仅在8、9日2 d可观察到明显的半日潮，之后受到多源内波的影响，2个站位的温度数据均出现无规律起伏，甚至覆盖掉半日潮的大周期。提取2个站位的25 ℃等温线及其包络如图3所示，从等温线包络线可以观察到，在试验初期2个站位的等温线变化趋势基本一致，之后出现较大差异，2个站位的25°等温线深度相差高达20 m，推测试验期间不仅存在B-O方向的内波，试验后期还在O站位出现了垂直等深线方向传播的内波，由于未观测到更大空间范围的温度变化情况，对于内波的运动方向和速度尚不能确定。

图2 9月8日-9月15日连续水文监测数据Fig.2 Thermistor string measured temperature from Sep 8th to Sep 15th

图3 B站位和O站位25 ℃等温线分布Fig.3 The 25 ℃ temperature contour of station B and station O

试验中采用CTD测量得到了全海深的温盐深数据，如图4所示。目前已经建立了许多利用温度、盐度和深度来计算声速的经验公式，为方便，采用Mackenzie[14]给出的9项算法，利用现场CTD测量得到的盐度数据代入经验公式，可将图2中的温度数据转化为声速数据，二者变化趋势基本一致：

图4 CTD测量得到的温度和盐度随深度的变化Fig.4 The temperature and salinity data measured by CTD

c=1 448.96+4.591T-5.304×10-2T2+

2.374×10-4T2+1.340(S-35)+1.630×

10-2D+1.675×10-7D2-1.025×

10-2T(S-35)-7.139×10-13TD3

(1)

1.2 浅水声速剖面的经验正交函数表示

利用EOF方法表示声速剖面可以大大降低声速剖面预报所需要的参数个数，基函数确定后，用有限的几阶系数即可进行重构。

根据一段时间内测量到的N条SSP数据，可通过求平均值的办法得到背景声速剖面：

(2)

则声速起伏可表示为：

(3)

设每条剖面离散到M个深度上，则N个Δc1(z)可以组成一个M×N的声速起伏矩阵：

(4)

求得C的协方差矩阵：

(5)

对矩阵R进行特征值分解：

RU=Uλ

(6)

(7)

式中：um为要提取的第m阶经验正交函数，将M个特征值λm降序排列，选取前d个特征值对应的特征向量用来表示声速剖面：

(8)

式中αm为第m阶EOF系数。

通常情况下，在掌握了背景声速剖面数据并提取到EOF基函数后，式(8)就可以较为准确地描述声速剖面。同理，也可以由声速剖面、EOF基函数和背景声速剖面数据计算得到对应的EOF系数。

1.3 试验数据的统计特征分析

由1.1节的试验数据初步分析可知，在试验最后一段时间监测到了内波活动，因此选取试验最后一天(9月15日)的声速数据进行EOF分析，截取时间段为15日5∶00-15∶00，该时间段的声速起伏变化如图5所示，其中白色分割线划分出了训练集以及待预测的测试集，整个时间段内的声速剖面集称为全集。图6为使用图5中训练集声速数据和全集声速数据得到的前2阶EOF基函数和背景声速剖面，二者形状基本一致。分别利用2个站位分解得到的EOF基函数和背景声速剖面计算得到O站位声速剖面的EOF系数随时间变化如图7所示，2种方法得到的系数变化趋势基本一致。从O站位剖面中任取一条进行重构，重构剖面与原始剖面吻合较好，说明在全集时间段内，声速剖面基函数未发生较大变化，可使用训练集提取的背景剖面和EOF基函数重构测试集声速剖面而不产生系统偏差。

图5 参与统计分析的声速数据Fig.5 Sound speed data used for statistical analysis

图6 分别利用训练集和全集声速数据得到的前两阶EOF基函数和背景声速Fig.6 The two first modes and the background sound speed profile extracted from training sets and complete sets respectively

对比图7中第1阶EOF系数与图5中的温跃层深度的变化规律，可以发现，二者变化趋势基本一致，在温跃层深度变化较大的时刻，第2阶EOF系数较大，尤其在15日13∶00左右最为明显。Gaelle曾在研究中给出了EOF系数的物理解释[15]，第1阶EOF系数α1表示了温跃层的垂直位移，即α1越大则代表温跃层越浅，跃层深度的变化周期也反映在α1的演化趋势中，第2阶EOF系数α2则表示了温度梯度的改变，在温跃层变化剧烈的时刻，α2取值也比较大。由此可推论，前两阶EOF系数之间并非完全独立，二者取值存在某种内在联系，即α2=f(α1)。

图7 训练集提取的EOF系数(9月15日)Fig.7 EOF coefficient extracted from the training set (sep.15)

将训练集和测试集的前2阶EOF系数分别绘制成散点图如图8所示，图中散点对应EOF系数分布，实线是通过多项式拟合得到的曲线，训练集和测试集数据均呈现出开口向上的抛物曲线分布，这与文献[15]中墨西拿海峡试验中得到的 “新月形”分布是一致的。图8中同时给出了曲线拟合的一元二次多项式，且开口大小基本一致。图9给出了2条拟合曲线的比较，在第1阶EOF系数取值[-20, 40]，二者基本吻合，只在[-40，20]有偏差，最大偏差为4，由于第2阶EOF基函数贡献相对较小，认为这个偏差在可以容忍的范围。对于没有明显内波存在的同时刻B站位声速数据(见图10)，EOF系数之间的统计特征则显得杂乱无章，无法得到较精准的拟合曲线。这也进一步印证了可以从一段时间内EOF系数之间的统计特征来监测内波的存在。另一方面，可以将内波存在时的特殊统计模型作为先验输入，从而简化声速剖面的预测。

图8 有内波存在时声速数据的统计特征Fig.8 Statistical characterization of SSP in the presence of interval waves

图9 训练集和测试集拟合曲线比较Fig.9 Comparison of fitting curves of training and test sets

2 基于LSTM的声速剖面预测

2.1 LSTM技术原理

循环神经网络(RNN)是以序列数据为输入，在序列演进方向进行递归且所有节点链式连接的递归神经网络。因其具有记忆性、参数共享等特点，在对序列的非线性特征学习建模方面具有一定的优势，被用于各类时间序列预报。然而RNN在处理时间序列上距离较远的节点时会涉及到雅可比矩阵的多次相乘，带来梯度消失或者梯度膨胀的问题。

LSTM是基于RNN 的一种完善，增加了记忆单元输入门、输出门和遗忘门等机制，网络结构示意图如图11所示。各记忆单元内部运算机制为：

图10 与图5同时刻的B站位声速数据的统计特征Fig.10 Statistical characterization of SSP of station B at the same time of Fig.5

图11 LSTM模型结构Fig.11 The schematic diagram of LSTM

(9)

式中：It、Ft和Ot分别代表输入门、遗忘门和输出门向量值；W1、W2和b分别表示输入节点与隐藏节点、隐藏节点与输出节点的连接权值以及对应的偏置项；xt对应t时刻的输入；ht-1对应t-1时刻的输出，代表了LSTM的隐藏状态；σ为sigmoid激活函数，可将实数映射到[0,1]，1表示上一时刻单元的信息全部保留，0表示上一时刻单元的信息全部遗忘。这些“门”的结构实质上是使用sigmoid神经网络与按位相乘结合，使用sigmoid作为激活函数的全连接神经网络层会输出一个0～1的数值描述当前通过该门的信息量，使得自循环的权重是变化的，避免了梯度消失等问题，可以控制不同时刻信息的传递来提升整个网络对于长序列数据的处理能力。

2.2 试验数据分析

考虑到内波统计特征，只将第1阶EOF系数作为预测值，同时从训练集提取的特征还包括EOF基函数、背景声速剖面及内波统计特征。训练集共包含7 h的声速数据，从中提取第1阶EOF系数并将数据预处理为具有零均值和单位方差的标准化数据，最终输入LSTM网络进行训练，最大迭代次数为200，学习率设置为0.005，在150次迭代以后学习率降低为0.000 1，训练在200代截至时已达到收敛，图12为去标准化后的网络输出与实际测量值的比较，预测值变化趋势与实际测量值基本一致，均方误差为3.4，只在跃层深度剧烈变化的时刻，预测误差略有增加。

图12 LSTM对第1阶EOF系数网络预测结果Fig.12 The first expansion coefficient predicted by LSTM model

利用内波统计特征和预测得到的第1阶EOF系数即可计算得到第2阶EOF系数的变化规律，如图13所示，与测量值基本符合。进一步利用先验EOF基函数和背景声速剖面重构得到测试集声速剖面如图14所示，预测值与测量值基本一致，最大绝对误差为3.69 m/s, 图15为时间序列预测的均方根误差(RMSE)，最大误差为1.3 m/s。

图13 基于声速剖面统计特征的LSTM网络时前两阶EOF系数预测结果Fig.13 The two first expansion coefficients predicted by LSTM model based on the statistical characterization of SSP