0 引　言

随着长续航力自主水下机器人、水下滑翔机、波浪滑翔机、无人帆船等面向海洋环境观测应用的长续航力海洋机器人（Marine Vehicle）技术逐步成熟^[1]，海洋机器人集群构建可控机动、组网协同、智能自主、持久在位的机器人化区域海洋环境观测网络成为研究和应用的重要发展方向。相对于海洋环境的时空尺度，海洋机器人观测网络获取的温度、盐度、海流流速等海洋环境观测数据通常较为稀疏。如何基于稀疏观测数据实现观测区域海洋环境温度、盐度、海流流速等海洋环境要素空间分布场的高分辨率预测，以满足军民等多方面应用对高分辨率海洋环境信息保障的需求，是海洋机器人观测网络应用需要研究和解决的一项重要问题。

针对该问题，国际上已开展诸多相关研究工作。在客观分析^[2]、克里金^[3]、高斯过程回归^[4]等基于最优插值和统计估计等方法实施基于观测数据预测未观测位置的数据之外，建立海洋环境预测模型并通过同化观测数据估计模型结构与参数、修正模型预测输出，进而通过模型实现海洋环境状态高分辨率预测是该研究中的另一种主要方法。如文献 [5]～文献[7] 分别基于径向基函数、本征正交分解、动力学模态分解方法建立海洋环境预测模型，利用海洋机器人观测数据估计模型参数，进而使模型逼近海洋环境状态并实现准确的海洋环境预测。除上述方法之外，随着深度学习技术的发展，利用具有更强非线性函数逼近能力的人工神经网络建立海洋环境预测模型、实现数据驱动的海洋环境预测成为近年来热点研究方向。SHI等 ^[8] 提出用于空间场时间序列预测的卷积长短期记忆网络ConvLSTM已成为诸多海洋环境空间场时间序列预测研究的重要工具。

目前在基于深度神经网络模型实现海洋环境预测的研究工作中，主要通过历史数据训练海洋环境空间场时间序列预测模型^[9]。本文在目前研究工作基础上，面向基于海洋机器人稀疏观测数据实现高分辨率、准确的海洋环境要素空间场预测需求，研究利用海洋机器人观测数据实现神经网络预测模型的在线学习。本文应用ConvLSTM建立高分辨率区域海洋环境预测模型，在应用历史数据训练的神经网络初始预测模型基础上，利用海洋机器人稀疏观测数据和序贯数据同化方法进一步优化模型参数，使模型不断逼近实时海洋环境状态、预测输出不断逼近海洋机器人观测。针对大尺度高空间分辨率预测带来的高计算复杂性问题，提出一种非均匀空间分辨率预测策略实施海洋环境预测，有效降低计算资源需求。利用中国科学院沈阳自动化研究所研发的“海翼”水下滑翔机，在南海北部实施的区域海洋环境观测实验中获得的观测数据，进行了区域海洋环境二维水平空间温度场预测实验，验证了方法的有效性。

1 海洋环境预测模型 1.1 空间场时间序列预测

本文研究二维水平空间海洋环境要素空间场时间序列预测。给定空间区域 $\mathcal{F}=\left\{{\mathit{s}}_{\mathcal{F}}|{\mathit{s}}_{\mathcal{F}}\in {R}^{2}\right\}$，针对温度、盐度、海流流速等海洋环境要素，实现其高分辨率空间场 $\mathit{x}\left(\mathcal{F}\right)=\left\{x\left({\mathit{s}}_{\mathcal{F}}\right)|{\mathit{s}}_{\mathcal{F}}\in \mathcal{F}\right\}$ 的时间演化预测，在海洋环境信息保障等领域具有重要需求。为此，本文研究如式 (1) 所示的海洋环境要素空间场时间序列预测问题，在给定预测时刻前一组多个离散时刻要素空间场的最优估计的条件下，计算在预测时刻要素空间场的预测结果：

${\mathit{x}}_{k+1}^{f}=\underset{{\mathit{x}}_{k+1}}{\mathit{arg}{\text{max}}}{p}\left({\mathit{x}}_{k+1}|{\mathit{x}}_{k}^{a},{\mathit{x}}_{k-1}^{a},\cdots ,{\mathit{x}}_{k-L+1}^{a}\right)。$

(1)

式中：${\mathit{x}}_{k+1}$ 与 $\left\{{\mathit{x}}_{k},{\mathit{x}}_{k-1},\cdots ,{\mathit{x}}_{k-L+1}\right\}$ 分别为在预测时刻 ${t}_{k+1}$ 及预测时刻前 $L$ 个离散时刻的温度、盐度、或海流流速等海洋环境要素空间场；${p}\left(\cdot \right)$ 为关于 ${\mathit{x}}_{k+1}$ 的概率密度函数；$\left\{{\mathit{x}}_{k}^{a},{\mathit{x}}_{k-1}^{a},\cdots ,{\mathit{x}}_{k-L+1}^{a}\right\}$ 为 $\left\{{\mathit{x}}_{k},{\mathit{x}}_{k-1},\cdots ,{\mathit{x}}_{k-L+1}\right\}$ 对应的最优估计；${\mathit{x}}_{k+1}^{f}$ 为 ${\mathit{x}}_{k+1}$ 的预测。

为实现选定海洋环境要素空间场的预测，将二维水平空间区域 $\mathcal{F}$ 离散化为 ${N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}$ 个预测网格节点 $\mathcal{D}=\left\{{\mathit{s}}_{\mathcal{D},i}|i=1,\cdots ,{N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}\right\}$。式中：${N}_{{\mathrm{lat}}}$ 与 ${N}_{{\mathrm{lon}}}$ 分别表示在预测区域范围内对空间区域离散化后的纬向和经向的预测网格节点数量，从而将对空间场 $\mathit{x}\left(\mathcal{F}\right)$ 的预测转化为对 ${N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}$ 个离散网格节点上的海洋环境要素值 $\mathit{x}\left(\mathcal{D}\right)=\left\{x\left({\mathit{s}}_{\mathcal{D}}\right)|{\mathit{s}}_{\mathcal{D}}\in \mathcal{D}\right\}$ 的预测。图1以海水温度为例展示了对空间区域的网格化离散结果。

本文通过数据驱动方法，利用人工神经网络建立上述海洋环境要素空间场预测模型。

1.2 ConvLSTM模型

文献 [8] 提出的卷积长短期记忆网络模型ConvLSTM是实现上述预测任务的有效技术途径。ConvLSTM融合了卷积神经网络（CNN）和长短期记忆网络（LSTM）的优点，通过在输入到状态和状态到状态转换的过程中引入卷积结构，对传统全连接LSTM进行了扩展，从而能够有效地捕获高维数据中的时空相关性。ConvLSTM网络的单元结构如图2所示。

ConvLSTM的关键方程如下^[8]：

${{i}}_{k} = {\mathrm{sigmoid}}\left({\boldsymbol{W}}_{xi}\mathcal{*}{{x}}_{k} + {\boldsymbol{W}}_{{H}{i}}\mathcal{*}{{H}}_{k-1} + {\boldsymbol{W}}_{{C}{i}}\circ {{C}}_{k-1} + {\boldsymbol{b}}_i\right)，$

(2)

${{f}}_{k} = {\mathrm{sigmoid}}\left( {\boldsymbol{W}}_{xf}\mathcal{*}{{x}}_{k} + {\boldsymbol{W}}_{{H}{f}}\mathcal{*}{{H}}_{k-1} + {\boldsymbol{W}}_{{C}{f}}\circ {{C}}_{k-1} + {\boldsymbol{b}}_{f} \right)，$

(3)

${{C}}_{k} = {f}_{k}\circ {{C}}_{k-1} + {{i}}_{k}\circ \text{tanh}\left({\boldsymbol{W}}_{xC}\mathcal{*}{{x}}_{k} + {\boldsymbol{W}}_{{HC}}\mathcal{*}{{H}}_{k-1} + {\boldsymbol{b}}_{C}\right)，$

(4)

${{o}}_{k} = {\mathrm{sigmoid}}\left({\boldsymbol{W}}_{xo}\mathcal{*}{{x}}_{k} + {\boldsymbol{W}}_{Ho}\mathcal{*}{{H}}_{k-1} + {\mathit{W}}_{\mathcal{C}\mathcal{o}}\circ {{C}}_{k}+{\boldsymbol{b}}_{o}\right)，$

(5)

${{H}}_{k}={{o}}_{k}\circ \text{tanh}\left({{C}}_{k}\right)。$

(6)

式中：${i}$、${f}$、${o}$ 分别为输入门、遗忘门和输出门；${H}$、${C}$ 分别为存储短期记忆的隐藏层状态和存储长期记忆的细胞单元状态；$\boldsymbol{W}$ 与 $\boldsymbol{b}$ 分别为神经网络的权值矩阵与偏置项；$\mathrm{*}$ 为卷积运算；$\circ$ 为哈达玛积运算；${\mathrm{sigmoid}}(\cdot)$ 为输入门、遗忘门与输出门选用的Sigmoid激活函数；$\text{tanh}(\cdot)$ 为隐藏层状态和细胞单元状态选用的双曲正切激活函数。

ConvLSTM适用于处理高维空间场的时序输入，以预测未来的高维空间场时序输出。

1.3 海洋环境预测模型设计

基于上述ConvLSTM单元模型，本文构建海洋环境要素空间场时间序列预测模型 $NN\left(\cdot \right)$：

${{x}}_{k+1}^{f}=NN\left({{x}}_{k}^{a},{{x}}_{k-1}^{a},\cdots ,{{x}}_{k-L+1}^{a},{\boldsymbol{W}}_{k+1}^{f},{\boldsymbol{b}}_{k+1}^{f}\right)。$

(7)

式中：${\boldsymbol{W}}_{k+1}$ 与 ${\boldsymbol{b}}_{k+1}$ 分别为ConvLSTM模型结构中在 ${t}_{k+1}$ 时刻的权值与偏置；${\boldsymbol{W}}_{k+1}^{f}$ 与 ${\boldsymbol{b}}_{k+1}^{f}$ 分别为 ${\boldsymbol{W}}_{k+1}$ 与 ${\boldsymbol{b}}_{k+1}$在 ${t}_{k+1}$ 时刻的预测。通过历史数据和海洋机器人原位观测数据学习及估计神经网络预测模型 $NN\left(\cdot \right)$ 的模型参数，即可实现海洋环境要素空间场的时间序列预测。

基于ConvLSTM单元模型，本文设计的海洋环境要素二维空间场时间序列预测的神经网络模型结构如图3所示。

在每个时间步，ConvLSTM层的输入为一个张量，其维度为 ${N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}$。ConvLSTM层的卷积核尺寸设置为 $3\times 3$；定义卷积核数量为 ${N}_{f}$；填充模式为“same padding”。

ConvLSTM层输出序列中的最后一个时间步的输出特征维度为 ${N}_{\rm lat}\times {N}_{\rm lon}\times {N}_{f}$，该特征作为卷积层的输入。卷积层的卷积核尺寸为 $3\times 3$；卷积核数量为 $1$；填充模式为“same padding”，激活函数选用线性激活函数。卷积层的输出维度为 ${N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}$，与海洋环境空间场预测输出维度相应。

1.4 数据归一化处理

本文对训练神经网络预测模型的样本数据应用Z-score方法进行归一化处理，将原始海洋环境要素数据映射为均值为 $0$、标准差为 $1$ 的分布：

${x}_{i,j}^{Z} = \left({x}_{i,j} - {\bar{x}}_{j}\right)/{\delta }_{j} ，i = 1,\cdots ,{N}_{sa}，j = 1,\cdots ,{N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}。$

(8)

式中：${N}_{sa}$ 为训练样本数量；${x}_{i,j}$ 为第 $i$ 个训练样本在 ${\mathit{s}}_{\mathcal{D},j}$ 位置的要素值；${\bar{x}}_{j}$ 为式(9) 所示训练样本在 ${\mathit{s}}_{\mathcal{D},j}$ 位置的样本均值；${\delta }_{j}$ 为 式(10) 所示训练样本在 ${\mathit{s}}_{\mathcal{D},j}$ 位置的样本标准差；${x}_{i,j}^{Z}$ 为 ${x}_{i,j}$ 的归一化值。

${\bar{x}}_{j}=\frac{1}{{N}_{{{sa}}}}\sum _{i=1}^{{N}_{{{sa}}}}{x}_{i,j}。$

(9)

${\delta }_{j}=\sqrt{\frac{1}{{N}_{{{sa}}}-1}\sum _{i=1}^{{N}_{{{sa}}}}{\left({x}_{i,j}-{\bar{x}}_{j}\right)}^{2}}。$

(10)

通过该方法，有助于消除数据间的尺度差异、避免模型受到数据尺度影响，使模型更容易学习数据的特征和模式，在训练过程中有助于加速收敛，提高模型的训练效率和预测性能。

在应用模型进行海洋环境要素空间场预测时，对模型输入的要素空间场时间序列首先进行如下归一化处理：

${x}_{k - i,j}^{{\mathrm{in}}} = \left( {x}_{k - i, j}^{a} - {\bar{x}} _{j} \right)/{\delta } _{j} ，i = 0, \cdots ,L - 1，j = 1, \cdots ,{N}_{{\mathrm{lat}}} \times {N}_{{\mathrm{lon}}}。$

(11)

式中：${x}_{k-i,j}^{a}$ 为 ${t}_{k-i}$ 时刻在 ${\mathit{s}}_{\mathcal{D},j}$ 位置要素值的最优估计；${x}_{k-i,j}^{{\mathrm{in}}}$ 为模型输入时间序列中第 $i$ 时间步的输入在 ${\mathit{s}}_{\mathcal{D},j}$ 位置的归一化值。对于模型输出数据，通过式 (12) 所示的基于上述归一化过程的逆过程，得到海洋环境要素预测值：

${x}_{k+1,j}^{f}={x}_{k+1,j}^{{\mathrm{out}}}\times {\delta }_{j}+{\bar{x}}_{j} \text{，} j=1,\cdots ,{N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}。$

(12)

式中：${x}_{k+1,j}^{{\mathrm{out}}}$ 为模型在 ${\mathit{s}}_{\mathcal{D},j}$ 位置的输出值；${x}_{k+1,j}^{f}$ 为对应的要素预测值。

2 基于海洋机器人观测数据的模型学习

由于海洋环境特性动态变化，为了应用前述建立的神经网络模型在海洋环境要素空间场连续预测过程中保持较高的预测准确性，还需要通过观测数据在线学习神经网络模型参数 $\boldsymbol{W}$ 与 $\boldsymbol{b}$。本文将该神经网络模型参数视为系统状态变量，通过基于集合卡尔曼滤波^[10]的序贯数据同化方法，利用海洋机器人实测数据学习神经网络模型参数。

2.1 状态方程

在神经网络模型中，权重 $\boldsymbol{W}$ 与偏置 $\boldsymbol{b}$ 决定模型的行为和输出，本文通过海洋机器人观测数据在线学习和调整 $\boldsymbol{W}$ 与 $\boldsymbol{b}$。定义预测模型状态向量：$\boldsymbol{\alpha }= \left[{w}_{1},\cdots , {w}_{{N}_{w}},{b}_{1}\cdots ,{b}_{{N}_{b}}\right]$，其中，$w$ 与 $b$ 分别表示 $\boldsymbol{W}$ 与 $\boldsymbol{b}$ 中的元素，${N}_{w}$ 与 ${N}_{b}$ 分别表示对应的元素数量。

为了实施滤波，本文将 $\boldsymbol{\alpha }$ 的时间演化建模为如下所示的噪声驱动过程：

${\boldsymbol{\alpha }}_{k}^{f}={\boldsymbol{\alpha }}_{k-1}^{a}+{\boldsymbol{\eta }}_{k}。$

(13)

式中：${\boldsymbol{\alpha }}_{k-1}^{a}$ 为 ${t}_{k-1}$ 时刻模型状态的最优估计；${\boldsymbol{\alpha }}_{k}^{f}$ 为 ${t}_{k}$ 时刻基于 ${\boldsymbol{\alpha }}_{k-1}^{a}$ 的模型状态预测；${\boldsymbol{\eta }}_{k}$ 为服从零均值高斯分布的过程噪声。

2.2 观测方程

建立观测方程，用于通过滤波在线学习模型状态。在 ${t}_{k}$ 时刻，海洋机器人观测网络在 ${N}_{\mathcal{S}}$ 处空间位置$\mathcal{S}=\left\{{\mathit{s}}_{\mathcal{S},i}|i=1,\cdots ,{N}_{\mathcal{S}}\right\}$ 获得海洋环境观测数据 $\left\{{y}_{i}|i=1,\cdots ,{N}_{\mathcal{S}}\right\}$，表达为观测向量：${\mathit{y}}_{k}={\left[{y}_{1},\cdots ,{y}_{{N}_{\mathcal{S}}}\right]}^{{\mathrm{T}}}$，其值为海洋环境要素在相应位置的实际值与海洋机器人传感器的测量噪声之和：

${\boldsymbol{y}}_{k}={\boldsymbol{x}}_{k}+{\boldsymbol{\varepsilon }}_{k}。$

(14)

式中：${\boldsymbol{x}}_{k}={\left[{x}_{k}\left({\mathit{s}}_{\mathcal{S},1}\right),\cdots ,{x}_{k}\left({\mathit{s}}_{\mathcal{S},{N}_{\mathcal{S}}}\right)\right]}^{\text T}$ 为 ${t}_{k}$ 时刻在全部观测位置的海洋环境要素实际值。${\mathit{\varepsilon }}_{k}、~{\mathcal{N}}_{{N}_{\mathcal{S}}}\left({\bf 0},{\sigma }_{\varepsilon }^{2}\boldsymbol{I}\right)$ 为 ${t}_{k}$ 时刻服从零均值高斯分布、协方差矩阵为 ${\sigma }_{\varepsilon }^{2}\boldsymbol{I}$ 的观测噪声，其中 ${\sigma }_{\varepsilon }^{2}$ 为观测噪声方差，$\boldsymbol{I}$ 为 ${N}_{\mathcal{S}}$ 阶单位矩阵。

在式 (14) 中，${\boldsymbol{x}}_{k}$ 可以表达为：

${\boldsymbol{x}}_{k}={h}_{k}\left({\boldsymbol{\alpha }}_{k}\right)={GPR|}_{\boldsymbol{S}}\left(NN\left({\boldsymbol{\alpha }}_{k}\right)\right)。$

(15)

式中：${h}_{k}\left(\cdot \right)$ 为 ${t}_{k}$ 时刻的非线性观测算子，能够在 ${N}_{\mathcal{S}}$ 个海洋机器人观测位置将海洋环境要素值从模型状态空间映射到观测空间。本文将其具体分解为神经网络模型预测输出和高斯过程回归估计2个过程：其中 $NN\left({\boldsymbol{\alpha }}_{k}\right)$ 为神经网络模型利用 ${t}_{k}$ 时刻的模型状态在预测节点上输出的预测值；${GPR|}_{\mathcal{S}}\left(\cdot \right)$ 利用模型在预测节点上的输出结果，通过高斯过程回归方法估计所有海洋机器人观测位置的要素值。

基于以上方程，可得到如下形式的观测方程：

${\boldsymbol{y}}_{k}={{h}}_{k}\left({\boldsymbol{\alpha }}_{k}\right)+{\boldsymbol{\varepsilon }}_{k}。$

(16)

2.3 基于集合卡尔曼滤波的参数学习

如式 (16) 所示，观测方程为非线性方程，因此针对该非线性系统模型，本文基于集合卡尔曼滤波方法实现模型参数学习。其具体过程如下：

共考虑 ${N}_{{{em}}}$ 个集合成员，对初始状态向量施加随机扰动生成状态向量的初始集合为$\left\{{\boldsymbol{\alpha }}_{0,i}\right.|\left.i=1,\cdots ,{N}_{{{em}}}\right\}$。

对于集合中各成员，有如下式的状态演化方程：

${\boldsymbol{\alpha }}_{k,i}^{f}={\boldsymbol{\alpha }}_{k-1,i}^{a}+{\boldsymbol{\eta }}_{k,i} \text{，} i=1,\cdots ,{N}_{{{em}}}。$

(17)

式中：$\left\{{\boldsymbol{\alpha }}_{k-1,i}^{a}|i=1,\cdots ,{N}_{{{em}}}\right\}$ 为 ${t}_{k-1}$ 时刻的状态向量后验估计集合；$\left\{{\boldsymbol{\alpha }}_{k,i}^{f}|i=1,\cdots ,{N}_{{{em}}}\right\}$ 为 ${t}_{k}$ 时刻的状态向量预测集合；$\left\{{\boldsymbol{\eta }}_{k,i}~{\mathcal{N}}_{{N}_{\alpha }}\left({\bf 0},\boldsymbol{\tau }\right)|i=1,\cdots ,{N}_{{{em}}}\right\}$ 为 ${t}_{k}$ 时刻服从零均值高斯分布、协方差矩阵为 $\boldsymbol{\tau }$ 的过程噪声集合。式中：${N}_{\alpha }={N}_{w}+{N}_{b}$ 为状态变量总数。

${t}_{k}$ 时刻状态预测的误差协方差矩阵为：

${\boldsymbol{P}}_{k}^{f}=\frac{1}{{N}_{{{em}}}-1}\sum _{i=1}^{{N}_{{{em}}}}\left({\boldsymbol{\alpha }}_{k,i}^{f}-{\bar{\boldsymbol{\alpha }}}_{k}^{f}\right){\left({\boldsymbol{\alpha }}_{k,i}^{f}-{\bar{\boldsymbol{\alpha }}}_{k}^{f}\right)}^{\text T}。$

(18)

式中，${\bar{\boldsymbol{\alpha }}}_{k}^{f}$ 为状态预测的集合平均：

${\bar{\boldsymbol{\alpha }}}_{k}^{f}=\frac{1}{{N}_{em}}\sum _{i=1}^{{N}_{em}}{\boldsymbol{\alpha }}_{k,i}^{f} 。$

(19)

对 ${t}_{k}$ 时刻的观测 ${\boldsymbol{y}}_{k}$ 施加随机扰动，生成 ${t}_{k}$ 时刻的观测集合 $\left\{{\boldsymbol{y}}_{k,i}|i=1,\cdots ,{N}_{em}\right\}$。对于状态集合中各成员，${t}_{k}$ 时刻的状态后验估计为：

${\boldsymbol{\alpha }}_{k,i}^{a}={\boldsymbol{\alpha }}_{k,i}^{f}+{\boldsymbol{K}}_{k}\left({\boldsymbol{y}}_{k,i}-{{h}}_{k}\left({\boldsymbol{\alpha }}_{k,i}^{f}\right)\right)。$

(20)

式中，${\boldsymbol{K}}_{k}$ 为 ${t}_{k}$ 时刻的卡尔曼增益矩阵：

${\boldsymbol{K}}_{k}={\boldsymbol{A}}_{k}{\left({\boldsymbol{B}}_{k}+{\sigma }_{\epsilon }^{2}\boldsymbol{I}\right)}^{-1}。$

(21)

式中：

${\boldsymbol{A}}_{k}=\frac{\displaystyle\sum _{i=1}^{{N}_{em}}\left({\boldsymbol{\alpha }}_{k,i}^{f}-{\bar{\boldsymbol{\alpha }}}_{k}^{f}\right){\left({{h}}_{k}\left({\boldsymbol{\alpha }}_{k,i}^{f}\right)-{{h}}_{k}\left({\bar{\boldsymbol{\alpha }}}_{k}^{f}\right)\right)}^{\text T}}{{N}_{em}-1}，$

(22)

${\boldsymbol{B}}_{k}=\frac{\displaystyle\sum _{i=1}^{{N}_{em}}\left({{h}}_{k}\left({\boldsymbol{\alpha }}_{k,i}^{f}\right)-{{h}}_{k}\left({\bar{\boldsymbol{\alpha }}}_{k}^{f}\right)\right){\left({{h}}_{k}\left({\boldsymbol{\alpha }}_{k,i}^{f}\right)-{{h}}_{k}\left({\bar{\boldsymbol{\alpha }}}_{k}^{f}\right)\right)}^{\text T}}{{N}_{em}-1}。$

(23)

${t}_{k}$ 时刻状态后验估计的误差协方差矩阵为：

${\boldsymbol{P}}_{k}^{a}=\frac{1}{{N}_{em}-1}\sum _{i=1}^{{N}_{em}}\left({\boldsymbol{\alpha }}_{k,i}^{a}-{\bar{\boldsymbol{\alpha }}}_{k}^{a}\right){\left({\boldsymbol{\alpha }}_{k,i}^{a}-{\bar{\boldsymbol{\alpha }}}_{k}^{a}\right)}^{\text T}。$

(24)

式中：${\bar{\boldsymbol{\alpha }}}_{k}^{a}$ 为状态后验估计的集合平均，有

${\bar{\boldsymbol{\alpha }}}_{k}^{a}=\frac{1}{{N}_{{{em}}}}\sum _{i=1}^{{N}_{{{em}}}}{\boldsymbol{\alpha }}_{k,i}^{a}。$

(25)

基于以上计算，${t}_{k}$ 时刻预测区域范围内海洋环境要素场的后验估计和相应的不确定性分别为：

${\boldsymbol{x}}_{k}^{a}={{h}}_{k}^{*}\left({\bar{\boldsymbol{\alpha }}}_{k}^{a}\right)={GPR|}_{\mathcal{D}}\left(NN\left({\bar{\boldsymbol{\alpha }}}_{k}^{a}\right)\right)，$

(26)

$COV\left({\boldsymbol{x}}_{k}^{a}|{\boldsymbol{y}}_{k}\right)={\boldsymbol{C}}_{k}+{\sigma }_{\epsilon }^{2}\boldsymbol{I}。$

(27)

式中：

${\boldsymbol{C}}_{k}=\frac{\displaystyle\sum _{i=1}^{{N}_{{{em}}}}\left({{h}}_{k}^{*}\left({\boldsymbol{\alpha }}_{k,i}^{a}\right)-{{h}}_{k}^{*}\left({\bar{\boldsymbol{\alpha }}}_{k}^{a}\right)\right){\left({{h}}_{k}^{*}\left({\boldsymbol{\alpha }}_{k,i}^{a}\right)-{{h}}_{k}^{*}\left({\bar{\boldsymbol{\alpha }}}_{k}^{a}\right)\right)}^{{\mathrm{T}}}}{{N}_{{{em}}}-1}。$

(28)

${h}_{k}^{*}\left(\cdot \right)$ 在 ${t}_{k}$ 时刻 ${N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}$ 个预测网格节点位置将海洋环境要素值从模型状态空间映射到 $\mathcal{D}$。

2.4 模型参数初始化

在上述参数学习过程中，$\mathit{\alpha }$ 的初始值对于后续学习具有重要应用。为了提高学习效率及参数收敛速度，本文应用历史数据和梯度下降方法训练神经网络，以获得较好的 $\mathit{\alpha }$ 初始值。

对于实现高分辨率海洋环境预测，通常缺乏相应高分辨率的海洋环境要素历史数据。针对该问题，本文应用数值海洋模式预报系统输出的历史数据以及高斯过程回归方法生成高分辨率数据用于神经网络训练，以获得 $\mathit{\alpha }$ 的初始值。

对于应用梯度下降方法训练神经网络，本文设置代价函数如下：

${\mathcal{J}}_{MSE}=\frac{1}{{N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}}\sum _{i=1}^{{N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}}{\left({x}_{pr,i}-{x}_{tr,i}\right)}^{2}。$

(29)

式中：${\mathit{x}}_{pr}$ 与 ${\mathit{x}}_{tr}$ 分别为模型预测输出结果和训练数据的实际值。本文应用Adam算法实现网络参数优化。

此外，对于海流流速空间场预测，考虑到海水的流动性并假设海水不可压缩，在预测区域范围内，一段时间内的海水净流入量应为 $0$。根据这一特点，海流流速空间场初始预测模型训练过程采用的代价函数中包含能够反映预测区域范围边界位置海流 $\left\{{x}_{pr,{e}_{j}}|j=1,\cdots ,2\times \left({N}_{\text{lat}}+{N}_{\text{lon}}-2\right)\right\}$ 流速净值的正则项，如下式：

$\begin{array}{c}{\mathcal{J}}_{MSE} = \displaystyle\frac{1}{{N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}}\displaystyle\sum _{i=1}^{{N}_{{\mathrm{lat}}}\times {N}_{{\mathrm{lon}}}}{\left({x}_{pr,i}-{x}_{tr,i}\right)}^{2} +\\ {\lambda} \left|\displaystyle\sum _{j=1}^{2\times \left({N}_{{\mathrm{lat}}} + {N}_{{\mathrm{lon}}} - 2\right)}{x}_{pr,{e}_{j}}\right|。\end{array}$

(30)

式中：${x}_{pr,{e}_{j}}$ 为第 $j$ 个边界位置海流值的模型输出结果；$\lambda > 0$ 为权重系数。

3 非均匀分辨率预测节点

目前海洋环境预测应用中，主要基于均匀空间分辨率网格节点划分实现预测。采用均匀空间分辨率网格节点划分，利用神经网络对大尺度海洋环境进行高分辨率预测，将需要计算大量节点进而导致庞大的网络规模和较高的计算资源需求，这对模型高效学习带来挑战。本文提出采用非均匀空间分辨率的预测网格节点划分策略，解决局部高分辨率预测与高计算资源需求之间的矛盾。

基于海洋环境要素空间分布的非均匀性，本文研究模型应用非均匀网格节点的环境预测，对关注区域进行加密预测，提高预测分辨率和准确性；对非关注区域进行稀疏预测，以实现在满足预测需求情况下降低模型复杂性和计算量，提升模型学习能力的目的。

为实施非均匀空间分辨率预测，首先需要确定分辨率调节目标，如对海洋环境不确定性较大或要素水平梯度变化剧烈区域实施加密预测。对于不同目标，本文生成非均匀空间分辨率预测网格节点的计算方法均相同，因此以对海洋环境要素随时间变化较大区域进行加密预测为例，给出如下所示非均匀节点位置的计算方法：

1）基于历史数据或神经网络模型预测输出数据，计算 $\mathcal{D}$ 内各节点位置方差，获得海洋环境要素方差分布场；

2）通过有限元求解过程中常用的自适应网格划分方法^[11]，利用三角面元覆盖方差场，对方差场进行三角网格剖分，在方差大的位置加密计算网格，并得到由 ${N}_{p}$ 个剖分网格顶点对应的点云分布 $\left\{{{p}}_{i}|i=1,\cdots ,{N}_{p}\right\}$；

3）采用聚类算法对点云进行聚类，获得设定预测节点数量 ${N}_{c}={N}_{\text{lat}}\times {N}_{\text{lon}}$ 的聚类中心 $\left\{{\boldsymbol{c}}_{j}|j=1,\cdots , {N}_{c}\right\}$，进而将所有聚类中心作为模型的预测节点位置。

本文采用模糊C-均值聚类算法^[12]实现聚类，通过最小化以下目标函数确定各聚类中心位置：

$\mathcal{J}\left({u},{c}\right)=\sum _{i=1}^{{N}_{p}}\sum _{j=1}^{{N}_{c}}{{u}}_{ij}^{2}{\left\|{{p}}_{i}-{{c}}_{j}\right\|}^{2}。$

(31)

式中：${{u}}_{ij}$ 为第 $i$ 个点云数据 ${{p}}_{i}$ 相对于第 $j$ 个聚类的隶属度；${\boldsymbol{c}}_{j}$ 为第 $j$ 个聚类的聚类中心位置；$\left\|\cdot\right\|$ 为欧氏距离。

4 实验验证

本节以区域海洋环境二维水平面温度场预测为例，应用中国科学院沈阳自动化研究所于2019年在南海北部布放的水下滑翔机观测数据，验证提出的区域高分辨率海洋环境预测方法的有效性。

4.1 预测区域

水下滑翔机作为一种浮力引擎驱动的海洋机器人，已成为执行长期大范围海洋环境观测任务的重要工具。

本文应用6台“海翼”水下滑翔机于2019年在南海北部执行区域海洋环境观测任务中获取的部分观测数据，进行本文提出方法的有效性验证。图4所示为观测实验区域位置及区域范围内水下滑翔机轨迹，轨迹旁所示编号如“1000J026”为该观测任务设置的水下滑翔机编号。根据如图4所示的实验区域，本文设置实施预测的区域海洋环境为北纬 20.80°～21.76° N 、东经119.76°～120.72° E 范围内 0.96° × 0.96°区域。在本实验中，应用该区域6台水下滑翔机在观测过程中通过温盐深传感器获取的100 m深度的温度数据进行该深度二维水平空间温度场预测。

4.2 参数设置

在该实验中，设置预测区域 $\mathcal{D}$ 内的预测网格节点数量为 ${N}_{c}={N}_{\rm lat}\times {N}_{\rm lon}=$11881，其中 ${N}_{\rm lat}={N}_{\rm lon}=109$。基于该预测网格节点数量，设置了相应的神经网络模型输入和输出张量维度。此外，神经网络预测模型输入的空间场时间序列长度设置为 $L=8$，ConvLSTM层的卷积核数量设置为 ${N}_{f}=8$。

为了确定预测节点的位置，本文基于HYCOM（Hybrid Coordinate Ocean Model）数值海洋模式预报输出数据计算预测区域范围内2019年05月14日08时前 24 h的方差分布场，并基于该分布场及本文采用的方法计算获得相应的非均匀预测网格节点，结果如图5所示。

可知，采用非均匀空间分辨率的预测网格节点，在温度场随时间变化较大的区域得到较高空间分辨率的节点，由此计算得到的最小值可达约为 0.036 km 。而在应用均匀网格条件下，若达到相同的空间分辨率，在计算过程中共需要计算约 ${N}_{c}\times {\left(111/ 112.5\right)}^{2}/ {0.036}^{2}\approx {N}_{c}\times 751$ 网格点。因此，通过非均匀空间分辨率的预测网格节点，可以在降低计算量的同时，达到在空间场随时间变化较大的区域实施高分辨率预测的目的。

为了提高模型学习效率，在同化水下滑翔机观测数据前，本文首先利用HYCOM数值海洋模式的历史输出数据训练神经网络实施模型初始化。HYCOM输出数据的均匀网格水平空间分辨率为0.08°，时间分辨率为 3 h 。为了实施如上所述的高分辨率区域海洋环境预测，本文通过高斯过程回归方法生成预测区域范围如图5所示的预测网格节点位置的高空间分辨率数据。在训练神经网络初始预测模型时，选用HYCOM数据训练样本的时间范围为2019年05月14日08时—08月12日05时，共对应 ${N}_{sa}=802$ 组训练样本。

在通过梯度下降方法训练得到的模型基础上，应用6台水下滑翔机的观测数据通过集合卡尔曼滤波方法学习模型参数。本实验过程共实施连续40次的学习，每次学习间隔时间为 3 h，每次学习时利用过去 3 h期间6台水下滑翔机获取的水下 100 m 深度的温度数据，该学习过程覆盖的时间范围为2019年08月12日08时—08月17日05时。学习过程中，模型状态变量数 ${N}_{\alpha }=$2697，状态变量集合成员数设置为 ${N}_{em}=100$，初始扰动协方差矩阵设置为 ${10}^{-5}\mathit{I}$，状态向量预测过程噪声协方差矩阵设置为 $\boldsymbol{\tau }={10}^{-5}\mathit{I}$，观测噪声方差设置为 ${\sigma }_{\varepsilon }^{2}={10}^{-4}$。

4.3 实验结果

为评估预测和学些结果，本文计算各学习周期在同化水下滑翔机观测数据后，基于 ${\boldsymbol{y}}_{k}={{h}}_{k}\left({\boldsymbol{\alpha }}_{k}\right)$ 计算的在水下滑翔机观测位置的预测数据，相对于水下滑翔机实际观测数据的均方根误差：

$RMS E\left({t}_{k}\right)=\sqrt{\frac{1}{{N}_{\mathcal{S}}}\sum _{i=1}^{{N}_{\mathcal{S}}}{\left({x}_{k}^{a}\left({\mathit{s}}_{\mathcal{S},i}\right)-{y}_{k,i}\right)}^{2}}。$

(32)

式中：${x}_{k}^{a}\left({\mathit{s}}_{\mathcal{S},i}\right)$ 为 ${t}_{k}$ 时刻预测模型在同化水下滑翔机观测数据后，通过式 (15) 计算得到的在水下滑翔机观测位置 ${\mathit{s}}_{\mathcal{S},i}$ 的预测结果；${y}_{k,i}$ 为对应水下滑翔机在 ${\mathit{s}}_{\mathcal{S},i}$ 的实际观测值。

由图6所示及式(32)计算得到的RMSE随学习周期的变化情况。结果表明，通过观测数据在线学习神经网络模型参数后，模型在水下滑翔机观测位置的预测结果相对于水下滑翔机观测的均方根误差由初始预测模型的 1.10℃，在10个学习周期后达到约 0.26℃，降低了76%。

图7为模型在2019年08月17日05时于预测范围内的海水温度场预测输出结果，图8为相应的预测方差场结果。

5 结　语

面向基于海洋机器人集群构建的区域海洋环境观测网络实施海洋环境信息保障需求，提出一种基于海洋机器人观测网络稀疏，观测数据预测高分辨率海洋环境要素空间分布场的方法。基于卷积长短期记忆神经网络，建立了有效的数据驱动区域海洋环境空间场时间序列预测模型，在基于历史数据训练的模型参数基础上，应用海洋机器人观测数据与集合卡尔曼滤波方法，实现了神经网络模型的在线学习，使其预测输出逐步逼近实际观测数据。针对大尺度海洋环境高分辨率预测所带来的计算复杂性问题，应用非均匀空间分辨率预测节点有效平衡了预测分辨率与计算量之间的矛盾。应用水下滑翔机网络实测数据实施的区域海洋环境二维水平面温度场预测实验，验证了提出方法的有效性。

在本文研究的方法中，海洋机器人在不同位置获取的观测数据对于模型学习效率、以及预测结果都具有重要影响。在目前研究结果基础上，优化海洋机器人的观测位置和路径，使其获取的观测数据能够最大化提高模型参数估计的准确性，是需要进一步开展的研究工作。