海冰工况下海洋运输系统的仿真建模分析

张 明

（福建船政交通职业学院，福建 福州 350007）

中国的渤海和北黄海海域地区是我国重要的经济开发区，其中渤海是我国的海上交通要道，是国内与国外沟通的海上门户之一，且目前已发现丰富的油气等资源，油气勘探、航海运输等系统正在快速的发展［1］。但是，每年的冬季，这些海域都会发生海冰，尤其是辽东海湾，海冰尤为严重。这些海冰的堆积和漂移对我国冬季的油气勘探尤其是海洋运输系统都造成了巨大的威胁。历史上曾经发生过海冰堵塞航道、撞毁海上运输船只以及推倒石油平台的事故［2］，为保证海洋运输系统的安全、适应冬季海洋油气的开发需求、减少并预防海冰可能会造成的灾害，亟需建立海冰工况综合预报系统。

目前，我国对于海冰的监测方式主要包括以下三种方式：一是日常的岸站监测，即海冰的实测数据以及气象数据进行预测。但是这种方式由于数据较为稀疏，往往以点窥面，误差较大；二是利用海洋、气象与海冰之间的统计学关系，按照以往的经验进行海冰的厚度、边界等的预测；三是利用渤海或黄海海域的热力学和动力学研究，预测海冰的增长、消融和漂移趋势［3］。这两种方式受制于数值计算本身的精度和效率，且具有一定的时效性，无法进行长期的预报。

近年来出现的遥测法，包括采用卫星遥感、无人机航测等方式搜集海冰相关资料，以及后来增加的船舶和野外观察也开拓了海冰的资料收集，这些可以作为海冰量值定标的依据。将这些数据与日常海冰监测方法结合，对海冰进行仿真建模分析，可以有效地为海洋运输系统提供依据。因此，本文将在海冰工况下对海洋运输系统的仿真建模进行分析。

1 海冰监测系统设计

1.1 海冰资料的收集

对海冰工况进行仿真建模分析，其中的一个重要环节是收集海冰相关资料。而海冰的研究与多个学科相关，包括海洋环境、大气等，所需要的资料包括卫星遥感资料、雷达监测资料、航测资料以及固定的平台资料等。该海冰工况仿真建模过程的设计简图如图1 所示。

图1 海冰工况仿真建模过程的设计简图

在进行海冰资料收集时，一般通过中分辨率成像光谱仪（MODIS）拍摄卫星遥感数据，以获取大面积的海冰信息。对遥感数据进行处理时，其反演流程如图2 所示。

图2 遥感图像的处理流程图

对小尺度范围的海冰图像或者需要海冰的细节信息，则采用在海冰监测站进行观察，监测信息主要包括会影响到海冰的气象、水文和海冰信息。气象信息包括风速、湿度、太阳辐射、气压和温度，气象信息每一秒测一次，可在监测站实时显示；水文信息包括水流速、流向、水文和盐分等，每隔15分钟采集一次，并在监测站显示；海冰信息包括海冰的类别、密集度、移动速度和移动轨迹等，这些采用雷达观测，海冰厚度采用CCD 摄像机进行观测。

2 海冰热力学和动力学分析

在对海冰进行数值计算之前，首先需要对决定海冰冻结和融化的热力学过程、决定海冰移动和变形的动力学过程、动力学和热力学共同决定的海冰厚度变化、质量和密集度的变化等进行研究［4-5］，才可对海冰进行后续的数值计算。

2.1 热力学过程

大气参数包括风速、相对湿度、气温、气压以及太阳辐射等，是海冰-海洋耦合的驱动条件。现在研究的海冰热力学模型都是在Stefan 公式的基础上进行研究得到，该公式认为海面结冰后，海冰表面的温度与气温近似，海冰底部的生长速度由结冰潜热、冰内部的热传导和冰内部的温度梯度决定，这种方法可以称为冰冻融冰度日法。但是，不同海域的海冰状况不同，但是基本原理相同，可以在Stefan 公式基础上改动得到渤海区域的冰厚计算公式，渤海区域的平整冰厚度hs［6-7］为：

式中α 为海冰的增长系数；FD 为每日平均气温低于-2℃的累积的海冰冰冻日期；TD 为每日平均气温高于0℃的累积的海冰解冻日期；KD为海冰刚出现时所需要的平均海冰冰冻日期。不同海域的海冰增长系数以及所需要的平均冰冻日期不同，在进行其他海域的海冰厚度计算时，将这两个系数更新即可。该公式忽略了太阳辐射以及海洋的热通量等对于海冰产生热影响的要素，因此在进行计算时，还需要将这些要素考虑进去。

在对海冰进行数值模拟仿真时，其热力学模型的基础是海冰内部的热力学平衡，这些热力学因素主要为冰冻融冰度日法中没有考虑到的因素，包括冰面上的太阳辐射、长波辐射、感热和潜热、海冰内部的热传导以及底部的海洋热能，海冰内部的热传导可以通过下式计算［8］：

式中ps，cs和ks分别为海冰的平均密度、比热容和热传导系数，Ts为海冰的平均温度；λ 和β 为根据经验确定的常数；Ss（z）为海冰内部距离海平面为z 深度的盐度；QSe为太阳到海冰上表面的有效辐射量；γs为海冰对太阳辐射的衰减系数；τs为海冰对太阳光的透射率；以上各公式中的下标s 表示海冰。

海冰的厚度变化在海冰的上下表面同时进行。对于海冰上表面来说，若海冰表面的温度大于海冰的融点，且表面热量的吸收和散发的总和大于0，则海冰表面融化；对于海冰的下表面来说，其变化主要由海冰的内部热传导、海洋热通量以及太阳辐射透射量来决定。将海冰的上表面和下表面同时考虑，可以得到海冰厚度的增长率为：

式中，下标u和下标d分别表示海冰上表面和下表面厚度的增长率，则海冰上表面和下表面厚度的增长率分别为：

式中hs为海冰的总厚度，与平均厚度不同，也可叫做物理冰厚；Ls为海冰融化时的潜热；Fw为海洋内部的热通量；Sd为到达海冰下表面的辐射透射量。

2.2 动力学过程

海冰在海洋上的漂移，其动力方程可通过牛顿定律表示，则单位面积海冰的动力学方程为：

式中M为海冰单位面积的质量，按照下式表示：

式中Ps和A分别为海冰的密度和密集度。为海冰在海上运动时的速度矢量；f为科氏参数，通过下式计算：

式中的ω为地球旋转的角速度；θ为海冰所处纬度位置。为与海面垂直的单位矢量；g为重力加速度；h0为海冰漂移处的海平面高度；为海冰内应力，当海冰密度较低或者与其他海冰无明显的相互作用，则该项可忽略。和分别为海冰处的风和水流对海冰的拖拽力，计算方式为：

式中Pf和Pl分别为海洋上空风以及海水的密度；Cf和Cl分别为风和海水对海冰的拖拽系数。

2.3 连续方程

海冰的质量变化是在热力学与和动力学共同作将用下引起的，海冰的连续方程可以将海冰的热力学和动力学过程耦合即可得到：

式中ΦA和Φh分别为海冰密集度、厚度的动力学变形系数；ΨA和Ψh分别为海冰密集度、厚度的热力学系数；vx和vy分别为海冰运动时在x和y方向的速度分量。

2.4 本构方程

在建立海冰的本构方程时，将海冰看做二维的连续体，并采用Hibler 粘塑性本构方程，该方程的建立是以符合椭圆屈服函数的二维粘塑性法则作为基础。该椭圆屈服函数为：

式中σ0和σ1为海冰内部的两个主应力；F0为海冰内部强度；e 为椭圆的主轴比。Hibler 方程中，海冰的应力和应变之间的关系如下：

式中σmn为海冰的应力张量；εmn为海冰的应变速率张量；δmn为Kronecker 算子；ζ为非线性块体的粘性系数；η为切向粘性系数。海冰内部的静水压力F0为：

式中F＇和A均为根据经验确定的常数。

3 海冰数值模拟计算

3.1 海冰数值模拟计算方法的选择

海冰在进行数值模拟计算时，一般将海冰看作是二维的流体或者流变体。针对流体的计算，按照采用坐标系的不同可以分为欧拉方法和拉格朗日方法。

欧拉方法适用于解决大变形的流体运动场。但对流场中包含两种以上的流体介质情况，介质之间的分界面以及混合网格中的各类运输量难以界定，使得人们不得不采取其他方法解决这类问题，海冰数值模拟即是这种情况。

拉格朗日方法适用于扭曲不严重的流体，尤其适用于多种介质的整体或局部运动。但是，若流场变形较大，计算过程中容易出现网格畸形和交叉的情况，使计算中断。对于海冰来说，其数值模型扰动较大且为多介质，可以将拉格朗日和欧拉方法结合进行计算，即网格质点的方法（PIC方法）。这种方法可以将网格的介质划分为若干质点，这些质点均有自己的质量、动量和能量，通过对质点进行运输即可实现对流体介质的运输。拉格朗日坐标下的光滑质点法（SPH）属于网格质点法（PIC 法），因此采用拉格朗日坐标下的光滑质点法（SPH）对海冰数值进行模拟计算。

3.2 海冰数值的模拟计算

采用拉格朗日坐标系的光滑质点法（SPH）对海冰进行数值模拟计算时，首先对海冰的动力方程进行处理可以得到：

其中海冰内力这项可以写为：

对于海冰中的第p个质点，具有以下性质：位置矢量rp，质量mj，质量密度Mp，海冰的平均厚度hp，质点处的应力张量σp，则海冰内力可以表示为：

将上式带入海冰的动力方程，并分别对动力方程的x和y向求解，可以得到p质点的动量平衡方程如下式所示：

4 神经网络对海冰面积的预报研究

海冰的数值计算可以实现海冰面积的短期预测，尽管预报结果具有一定的精度的，但仍然存在一定的误差。为能够及时、有效地识别海冰并进行海冰预警，以尽早采取措施降低海冰造成的损失，本文引入了预报精度更高的神经网络模型对海冰面积进行预报，同时为保证海冰面积在检测过程中的精度，采用卡尔曼滤波算法对预测结果进行跟踪，以进一步提高预报精度。

4.1 人工神经网络模型的设计

人工神经网络是一种数学模型，通过利用类似大脑神经突触的结构进行数据和信息的处理。在这些人工神经网络模型中，BP 神经网络或者其改进的形式是最常用的网络模型。根据本文的需求，适合对海冰的面积进行预测的模型只有前馈网络中的BP 神经网络。另一方面，BP 神经网络具有精度高、学习效率快、泛化能力强以及快速收敛的特点，因此本文采用BP 神经网络算法完成海冰面积的预测。

BP 神经网络模型在进行设计时，主要对其网络层数、每层的神经元个数、各层之间的激活函数、初始值和学习速率进行设计。对于BP 神经网络的输入值层和输出层，主要根据需要解决的问题以及最终的需求决定。根据设计要求，本文设计的BP 神经网络的输入层主要包括海冰冰冻温度、海冰融化温度、每日平均气温和离岸大风引起的温度变化，输出层为未来四年的预测海冰面积。

隐含层节点数量对BP 神经网络的预测结果影响较大，节点数量过多会使BP 神经网络过于复杂，延长网络训练以及后续计算的时间，缺少泛化能力，甚至降低识预测的准确率；而节点数量过少则会无法准确提取训练样本特征，无法准确识别未进行训练的样本，容错性差。在设计过程中，确定隐含层节点数量k 主要通过经验公式计算，如下式所示：

式中m和n分别为输入层和输出层的神经元个数，a为（1，10）范围的常数。

BP 神经网络在进行训练以及最终的结果预测之前，需要初始化神经网络参数，即模型的所有神经元之间的激活函数均设定小于1 的任意值，一般取（-1，1）之间的任意值。

其后选择BP 神经网络的训练样本，样本尽量选择输入层所包含的范围，才能获得良好的网络性能。根据以上原则，基本可以确定BP 神经网络的预测模型，训练步骤如图3 所示。

图3 BP 神经网络的预测步骤

图中的正、负样本分别指图像中包含和未包含海冰的图像。经过参数的调整和训练后，可以确定BP 神经网络的激活函数采用线性函数purelin，隐含层神经元个数为4，目标误差为0.01，最大的训练次数为1000。

4.2 卡尔曼滤波算法设计

在使用航测资料进行海冰面积的预测时，需要采用卡尔曼滤波算法对图像进行跟踪处理。算法对图像的求解过程如图4 所示。

图4 卡尔曼滤波算法对图像的求解过程

图像在t 时刻的先验估计值Xt|t-1根据滤波方程计算：

式中Zt为描述t-1 和t 时刻图像状态关系的矩阵；Kt为控制图像输入的矩阵；μt为控制系统输入的矩阵。其后根据先验估计值，计算该值对应的协方差矩阵，计算方式如下：

第三步为计算系统当前的卡尔曼增益Kt，计算方式如下：

式中Qt为描述t 时刻图像状态映射到无人机观测区域的矩阵；Gt为按照高斯分布表示的协方差矩阵。第四步为计算系统更新后的后验估计值Xt|t和协方差矩阵，分别如下式所示：

式中Zt为当前时刻观察的图像值。其后，系统将此时刻的图像信息作为上一时刻，依次循环该求解过程。

4.3 海冰的结果预测和分析

为检验模型的对海冰面积预测的精度，采用该模型对渤海海湾的海冰面积进行预测。对2019年12 月1 日至2020 年3 月8 日的海冰面积进行预测，同时采用遥感以及固定平台对海冰的实际面积进行观察，将预测值与实际的测量值进行对比，计算预测误差。结果如表1 所示。

表1 海冰预测试验结果对比

由表1 可知，对海冰面积的预测值和观察值的结果以及海冰面积的变化趋势基本一致，说明该模型可用于对海冰面积的预测，可以采用该模型为海洋运输系统的正常运行提供理论依据。

5 结论

当我国的渤海和北黄海海域发生海冰时，会对航海运输和油气勘探等造成破坏，为保证海洋运输系统的安全、减少并预防海冰可能会造成的灾害，本文主要进行了以下工作：

（1）建立了多个收集海冰相关资料的方法，包括卫星遥感、雷达监测、航测以及固定的监测站等。

（2）对海冰的热力学和动力学要素进行了分析。通过对冰面太阳辐射、长波辐射、感热、冰内热传导和冰下海洋热通量等热力学要素对海冰生消机理的影响进行分析，讨论海冰厚度存在的条件；通过对拖曳力、海冰内力等动力学因素进行分析，讨论了海冰漂移的动力学方程。同时建立了海冰的连续方程和本构方程。

（3）充分考虑了渤海海域海冰生消和漂移的特点，采用拉格朗日坐标系的光滑质点法对海冰的演化过程进行了数值模拟。

（4）本文采用BP 神经网络对海冰面积进行预报，并采用卡尔曼滤波算法对预测结果进行跟踪，以进一步提高预报精度。采用该模型对渤海海冰面积进行预测，结果表明该模型的精度较高，可用于对海冰面积进行预测，为海洋运输系统的正常运行提供理论依据。