进港航道中船舶搁浅概率的蒙特卡洛仿真

任亚磊，牟军敏，李亚军，易 侃

（1内河航运技术湖北省重点实验室（武汉理工大学），武汉 430063；2中华人民共和国上海海事局，上海 200086 3高性能船舶技术教育部重点实验室（武汉理工大学），武汉 430063）

进港航道中船舶搁浅概率的蒙特卡洛仿真

任亚磊1，2，牟军敏3，李亚军1，易 侃1

（1内河航运技术湖北省重点实验室（武汉理工大学），武汉 430063；2中华人民共和国上海海事局，上海 200086 3高性能船舶技术教育部重点实验室（武汉理工大学），武汉 430063）

文章通过对影响船舶吃水的主要因素，如船舶对波浪的响应、横倾引起的吃水变化、浅水中船舶的下坐以及航道底部等的分布规律的研究，拟合出其概率密度函数，然后使用蒙特卡洛（Monte Carlo）仿真船舶搁浅概率，其中船舶对波浪的响应采用了先进的粘性流数值计算方法。该模型可为受限水域航道建设和船舶安全通航提供参考。

蒙特卡洛；动态吃水；船舶搁浅概率

1 引 言

随着船舶的大型化和高速化，船舶搁浅的问题日益凸显。搁浅事故常带来船体结构损坏、环境污染甚至人员伤亡等严重后果，自上世纪50年代末以来国内外学者对其进行了一系列研究。祁恩荣[1]对前人的成果进行了总结，这些研究成果大概可归纳为三类：第一类主要从船舶实际操作入手开展研究，如采用实船事故调查、事件树、模糊事件树、决策树、贝叶斯网络等方法，探索搁浅的原因和机理，对船

舶偏离航道、机械失效、救援失效、恶劣海况等导致搁浅的因素进行分析，定性或定量地描述船舶搁浅风险的大小[2-5]；第二类主要从船体结构力学入手研究，采用缩尺或实尺度模型试验、有限元分析、简化数值模拟和解析经验公式等方法，研究搁浅事故中船舶力学问题，分析船底结构损伤、变形受力，为改善船体结构提供参考[6-7]；第三类则侧重于航道中船舶通航条件等相关数据的采集和分析，以充分利用航道水深，减小船舶搁浅风险。黄蕴和[8]针对重载船舶在进出浙江海门港乘潮作业问题，结合港外航道潮汐变化规律，研讨了现有航道水深条件下乘潮作业的合理性。胡勤友[9]综合潮汐数据、电子海图、船位和航向信息等开发了船用搁浅预警系统。

搁浅是个非常复杂的科学与工程技术问题，但其发生的必要条件比较简单，即航行水域水深相对船舶吃水不足。如大型船舶驶离设计航道导致搁浅，或在动静力作用下船舶浮态发生变化，导致在航道中吃水不足搁浅。考虑到上述文献中对动态吃水引起可航水深不足的研究相对偏少，本文主要围绕船舶动态吃水开展船舶搁浅概率的研究。总体上，影响船舶动态吃水的主要因素包括：潮汐、船舶对波浪的响应、浅水中船舶下坐、船舶操纵运动所致浮态变化以及航道底部情况等。

这些因素具有很强的随机性和不确定性。结合船舶水动力计算与航道工程技术，综合考虑以上各因素，已是目前搁浅概率研究的新方法。澳大利亚的O’Brien[10]博士创造了动态富裕水深系统（Dynamic Under-keel Clearance System）并开发了系列的产品，注册了DUKC的商标和专利。相比较传统的港口航道规划设计方法和建设使用过程中地毯式疏浚技术，这种应用系统融合了船舶波浪中运动响应、概率设计、气象预报及海洋水文等多要素信息，可以计算出一定置信度下船舶靠离港的最小安全富裕水深值和乘潮时间窗，充分利用了水深，减少了港口疏浚费用。另外一种较新的方法是蒙特卡洛（Monte Carlo）仿真。它在考虑影响船舶吃水各因素的基础上，根据其分布规律产生大量随机数，最后统计得到搁浅的频率，进而得出概率。Schoenmakers[11]采用该方法对开普敦港（Cape Town）进港航道的设计深度进行了验算，波浪处理采用了swan模型，船舶对波浪的响应采用了简化的经验公式，而且没有考虑船舶在波浪中的横倾和纵倾问题。Gucma[12]采用了类似的方法，更偏向于经验型，相关数据的获取采用了模拟器与实船实验相结合的方法，另外又考虑了可航的泥底、各种测量误差等因素，求出搁浅的概率。但受限于试验次数和实验条件的设置，误差等参数选取进行了大量简化。为此，本文根据蒙特卡洛仿真的思想，结合CFD计算结果，开展了设计代表船舶在航道中搁浅的概率计算，并以浙江省台州某港区航道进行了实例研究。

2 基本原理和方法

图1 导致船舶吃水变化的因素Fig.1 Causations to ship dynamic draft

蒙特卡洛方法，又称计算机随机模拟方法，是一种基于“随机数”的计算方法。该方法的基本思想就是用事件发生的“频率”来决定事件的“概率”。上世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得人们可以在计算机上利用数学方法大量、快速地产生一定统计特征的随机数进行模拟试验。目前这一方法已经广泛地运用到数学、物理、管理、生物遗传、社会科学等领域，并显示出了特殊的优越性。

实际应用中的随机数通常都是某些数学公式计算而产生的伪随机数，这些数不是严格意义上的随机数，但是，只要它们能够通过随机数的一系列统计检验，就可以将其作为真随机数使用。因此，理论上要求伪随机数产生器具备以下特征：良好的统计分布特性，高效率的伪随机数产生，伪随机数产生的循环周期长，伪随机数可以重复产生等。这些伪随机数产生器借助MATLAB等数学工具软件可以方便实现。船舶在航道中运动，引起吃水变化是高度复杂的过程，涉及了多种因素作用，如图1所示。

船舶所需吃水与航道水深的关系可以用（1）式表示：

其中：d 为航道深度（m）；Htide为潮高（m）；D 为设计船型吃水（m）；r为船舶对波浪的垂向响应（m）；Bf为航道底部情况（m）；h 为横倾引起的吃水变化（m）；s为浅水引起的下坐（squat）（m）。

从随机过程理论看，上式中许多因子是随机的，而且符合某种统计分布特征。因此根据蒙特卡洛方法求取船舶搁浅概率，即P（Z<0），可以通过研究上式中各因素的概率密度函数，产生大量的随机数序列，最后统计船舶搁浅的频率以代替搁浅概率。

值得一提的是，对于同一航道，不同的航段航行条件也有较大差异，如防波堤内外航道的风浪要素、以及航段曲率等。因此，为方便研究可根据不同的航行条件将航道划分为数段，分别进行研究。

3 相关参数求取

3.1 航道深度

对航道设计来说，航道深度是最终求取的目标。通常按照港口航道设计要求的搁浅率为临界点，或者遵循费效平衡的原则确定。

当然对于一定设计深度的航道，本文方法则可以计算出代表船型的搁浅概率，为船舶安全通航提供参考。

3.2 潮汐

潮汐因地而异，较复杂，但有其规律，可进行准确预报。合理利用潮汐，根据潮汐预报制作出船舶安全通过航道的时间表，可有效提高航道通过能力，减少疏浚工作量，节约成本。具体潮高可通过潮汐表求得。如果精确计算，应考虑统计船舶通过过程中潮汐变化；在实际操作中，若船舶通过某一航段所需时间较短，潮高变化不大，也可简化为采用此时间段内的平均潮高。

3.3 设计船型吃水

船舶吃水一般是指船舶的吃水深度，即船舶的底部至船体与水面相连处的垂直距离。船舶吃水有设计吃水和结构吃水两种，设计吃水是指船舶装载设计载重量货物情况下达到的吃水深度；结构吃水是船舶装载最大载重量货物情况下达到的吃水深度。这里取设计吃水作为检验值。

3.4 船舶对波浪的响应

船舶对波浪的响应是对波浪条件的统计响应，这里既有波浪要素的统计规律，又有船舶响应的统计规律。本文所述的船舶对波浪响应是对某种特定波浪条件的响应，即根据波浪要素的统计资料，采用设计波法描述工程水域的不规则波浪，亦即选取工程水域的不规则波特征要素构造规则波，以特定波浪条件来代替实际波，简化计算。

船舶在波浪中的运动与海域水深、潮位、波浪要素及船体自身的尺度和型线等直接相关。目前，实际海域中海洋结构物的波浪载荷及运动的研究方法主要有设计波法、设计谱法和时域计算方法等。将来随着仿真技术的进步，使用模拟器与实船实验相结合的方法获取数据将成为主流，不过现阶段受限于试验次数和实验条件，各参数选取都进行了大量简化，结果较粗糙，因此最终还是选用了CFD方法计算构造波中的船舶响应。

由于Boussinesq方程计算模型通常是将三维问题简化为平面二维问题进行计算[13]，本文采用粘性流数值计算方法，用体积分数VOF法跟踪自由面波形，计算船体在波浪中的粘性流场，得到船体受到的波浪力和船体的在波浪中的运动响应。计算船舶对波浪的响应的具体方法是[14]：

（1）采用设计波法描述工程水域的不规则波浪，即选取工程水域的不规则波特征要素，构造规则波。如选取若干年一遇为特征波重现期，分别选取其中十分之一波高H1/10和有义波高H1/3及其周期作为设计波的波高和周期，构造三阶Stokes规则波。

（2）选取代表船型船体周围水域为计算流体域，考虑水的粘性影响，通过边界造波法，求解粘性非定常Navier-Stokes（N-S）方程，再计算流体域生成三阶Stokes波，计算波浪中船体周围流场和船体受到的水动力。然后根据船体运动方程，得到船体在波浪中的六自由度运动。

3.5 航道底部情况

由于冲刷、淤积和施工工艺的不确定因素，航道底部不均匀等高。一般可以通过实际测量，拟合出其分布函数。

3.6 横倾引起的吃水变化

横倾是指船舶自正浮状态向左舷或右舷方向倾斜的一种浮态，此处仅考虑由于船舶转向操作引起的横倾，波浪引起的横倾的计算，在3.4节中已有介绍。船舶横倾导致吃水变化如图2所示，图中：C为船舶重心；α为横倾角；h为吃水变化。由此不难推出h的计算公式：

3.7 浅水引起的下坐（Squat）

船体在浅水中航行下坐，主要包括船体整体下沉和纵倾变化两部分。1967年Tuck利用细长体理论首先给出了一种估算方法。之后，在此基础上出现了许多解析式或半经验的计算公式[15]，如：Hooft公 式 （1974）、Huuska 公 式 （1976）、Eryuzlu 公 式（1978）、Barrass 公式 （1981）、Romisch 公式（1989）、Millward公式 （1990）、Millward 公式（1992）、Eryuzlu公式（1994）和 Ankudinov 公式（1996）。 但实际上，即使对同一种船型，用不同的公式计算差别仍然会较大，特别是在船速大于8节后，差别更为显著。因此，理论上任何一种船舶航行下沉量的计算方法都需要大量的实船测量数据进行验证和修正。如Gucma[12]采用GPS-RTK实验给多种模型赋权重值确定下坐。当然，在实测资料相对匮乏的情况下，实际计算可以根据各公式的应用范围、船型特点、航道特点以及通航密度等条件进行合理选取。而且从安全角度考虑，不妨采用下沉量较大值者。

一般来说航道中船舶的下坐可以采用BarassⅡ公式

图2 船舶横倾导致吃水变化示意图Fig.2 Ship heel and draft change

其中：Smax为由于浅水下坐导致的吃水变化（m）；Cb为船舶方形系数。

上式中：AS-vertical为船体横切面面积（m2），Ach为航道横切面面积（m2），VK为船舶的对地航速（kn）。

4 实例研究

4.1 工程概况

台州湾位于浙江中部沿海，为开敞的河口湾。海岸属于淤泥质或人工海岸，以平直的淤涨型岸滩为主。海湾面积为911.561 km2，其中潮滩面积258.748 km2，浅水区面积比例相对较高。随着海洋产业的发展，特别是世界级炼油和化工一体化工厂选址台州湾黄礁作业区，作业区港口航道远景规划5万吨级双向乘潮通航，全长12 km，底宽260 m，设计底标高-11.5 m。

4.2 数据处理

为方便研究根据航速、航道曲率、波浪大小等因素将航道划分为九段分别进行研究，如图3所示。

图3 航段划分示意图Fig.3 Sections divided for the approach channel

4.2.1 航道深度

航道远景规划5万吨级双向乘潮通航，全长12 km，底宽260 m，设计底标高-11.5 m。

4.2.2 潮汐

本海域无长期实测验潮资料，国家海洋局上海东海海洋勘测设计研究院在其南部松门镇横门设立潮位站，进行了一年的潮位观测（2008.10～2009.9）。根据实测验潮资料统计、计算得出本海域潮位特征值见表1。

表1 横门站潮位特征值Tab.1 Characteristics of tide level(Heng-men hydrographic station)

船舶通过航道时的具体潮高及其分布变化可通过潮汐表求得，本文是求取设计航道高水位时5万吨级油船搁浅概率，因此选取工程设计乘潮高水位5.33 m为典型水位进行简化计算，为提高仿真精度，应考虑潮高大于5.33 m时的潮汐水位变化分布，作为一个变量加到设计乘潮高水位上。

4.2.3设计船型吃水

所研究航道内行驶的主要运输船型是以1～3万吨级的成品油船为主，最大船型可达5万吨级。模型所选取标准船型主尺度如表2所示。

表2 标准船型主尺度Tab.2 Pariticulars of the designed ship

4.2.4 船舶对波浪的响应

本文选取的是波向频率最高的NE～ENE，以2年一遇为特征波重现期，以其H1/10作为设计波的波高和周期，构造三阶Stokes规则波。迎浪和横浪是船舶航行的两个极端条件，尤其是横浪对船体纵向和横向的浮态影响最大，本文计算了船舶迎浪和横浪下的纵摇、垂荡和横摇运动。波浪要素选取H1/10，即H=3.6 m，T=10 s，水深为h=16.83 m，船舶航速V=5 kns。计算结果如图4、5所示。

图4 迎浪船体垂荡、横摇及纵摇运动时历Fig.4 Ship response to head sea in time domain(heave,roll and pitch)

图5 横浪船体垂荡、横摇及纵摇运动时历Fig.5 Ship response to transverse sea in time domain(heave,roll and pitch)

4.2.5 航道底部情况

航道底部情况由于施工精度、冲刷、淤积等随机、不确定因素的影响，精确地描述非常困难。对于台州航道，通过测量，采用均值为0.1，标准差为0.1的正态分布近似描述。

4.2.6 横倾引起的吃水变化

此处只考虑船舶转向操作引起的横倾，该值大小与船型、船速、所操舵角、重心高度、初稳心高度等因素密切相关，实际观测表明其最大值一般不超过5°。Schoenmakers[11]认为船舶在弯曲航段，因操舵引起的横倾值一般不超过5°，具体的处理方法是通过航行仿真，选取每一航段的最大横倾值，数据选取从2.5°到5°不等；Gucma[12]则选取+/-3°简化处理。考虑到文献[11]中研究船型为集装箱船、文献[12]中研究船型为客船，而本文的研究对象为油轮，稳性大于集装箱船和客船，又结合几位船长的意见，最终选取操纵横倾为2°。

4.2.7 浅水引起的下坐（Squat）

根据BarassⅡ公式，求取下坐值关键在于获取船舶的方型系数和船速。所选取标准船型的方型系数参照表2，船速由各航段实测得出，取其平均速度，并以一定的分布规律进行描述，如航段1平均船速符合威布尔分布η=2.04，β=1.93；航段9符合均值为6.03，标准差为1.81的正态分布。

4.3 计算结果

分别以航段1（防波堤外直航道）和航段9（防波堤内弯曲航道）为代表进行计算，结果分别如图6至9所示。

从仿真结果看，10万次的仿真中，5万吨油船在航段1横浪条件下搁浅次数为49次，迎浪条件下为1次，航段9横浪和迎浪条件下都没有发生搁浅。结果表明，航道设计乘潮高水位16.83 m，具有较高的可靠性，可以满足设计船型在一般条件下的安全航行，且具有较大的安全裕量。如果以港口管理部门所能接受的搁浅率0.5%为临界点，则航段1的安全裕量为0.47 m，即航道设计深度为11.03 m即能满足设计要求。以5万吨油船每厘米吃水吨数65吨计算，现有航道水深可提升运力3 055吨。航段9的安全裕量为1.29 m，即航道设计深度为10.21 m即能满足设计要求。以5万吨油船每厘米吃水吨数65吨计算，现有航道水深可提升运力8 385吨。就航道疏浚而言，上述方法节省的疏浚量非常可观，经济效益明显。

图6 航段1横浪5万吨级满载油船搁浅概率仿真Fig.6 Probability Simulating for a 50 000 ton class laden oil tanker grounding in Section 1(transverse sea)

图7 航段1迎浪5万吨级满载油船搁浅概率仿真Fig.Probability Simulating for a 50 000 ton class laden oil tanker grounding in Section 1(head sea)

图8 航段9横浪5万吨级满载油船搁浅概率仿真Fig.8 Probability Simulating for a 50 000 ton class laden oil tanker grounding in Section 9(transverse sea)

图9 航段9迎浪5万吨级满载油船搁浅概率仿真Fig.9 Probability Simulating for a 50 000 ton class laden oil tanker grounding in Section 9(head sea)

5 结 论

从随机过程的角度考虑影响船舶吃水的主要因素：船舶对波浪的响应、横倾引起的吃水变化、浅水中船舶的下坐以及航道底部情况的分布规律，开展船舶搁浅概率的蒙特卡洛仿真，相对于传统基于规范的设计和分析方法，有其明显的优点。

本文在综合考虑水动力计算和航道工程技术的基础上，针对浙江台州某港区的设计航道水深进行了蒙特卡洛实例仿真，结果表明原规划航道设计水深还有较大的裕度。

[1]祁恩荣,崔维成.船舶碰撞和搁浅研究综述[J].船舶力学,2001,5(4):67-80.Qi Enrong,Cui Weicheng.A state-of-the-art review on ship collision and grounding[J].Journal of Ship Mechanics,2001,5(4):67-80.

[2]陈 刚,张圣坤.船舶搁浅概率的模糊事件树分析[J].上海交通大学学报,2002,36(1):112-116.

[3]于 鹏.搁浅船舶危险度定量研究[D].大连:大连海事大学,2011.

[4]Zhu Renqing,Zhao Hongjiang.Probability of ship grounding and collision calculation[J].Shipbuilding Science and Technology,1999(1):4-9.

[5]胡中凯,尹 群,刘海燕等.基于贝叶斯网络方法对船舶搁浅概率的研究[J].舰船科学技术,2010,32(2):23-26.

[6]刘 峰,王自力,崔维成.船舶结构的搁浅数值仿真研究[J].船舶,2006,3:24-27.

[7]杨传武,王安稳.小水线面船搁浅过程的数值仿真研究[J].固体力学学报,2006,29:24-27.

[8]黄蕴和,金丕信.海门港外航道通航条件研究[J].中国航海,1996,1:77-83.

[9]龚安祥,胡勤友,徐 铁.SAGA:一种新的船舶搁浅预警模型[J].大连海事大学学报,2007,28(1):106-110.

[10]O’Brien T.Experience using an innovative under keel clearance prediction system in australian ports[J].Port Tech.Int.,1999(6):165-169.

[11]Schoenmakers N W A,WIT A B F de.Probabilistic design entrance channel port of cape town[R].Minor thesis.Hydraulic Engineering,2006.

[12]Gucma L,Marta S.Monte Carlo method of ship’s under-keel clearance evaluation for safety of ferry approaching to Ystad port determination[J].Journal of Konbin,2008,5(8):36-44.

[13]王大国,邹志利,唐春安.港口非线性波浪耦合计算模型研究[J].力学学报,2007(5):587-594.

[14]武汉理工大学航运学院.台州湾岸线的浅水深用通航安全关键技术研究报告[R].浙江交通厅科技项目编号:2009W11,2009.

[15]洪碧光,于 洋.船舶在浅水中航行下沉量的计算方法[J].大连海事大学学报,2003,29(2):1-5.

Monte Carlo simulation for the grounding probability of ship maneuvering in approach channels

REN Ya-lei1，2,MOU Jun-min3,LI Ya-jun1,YI Kan1
(1 Hubei Key Laboratory of Inland Shipping Technology(Wuhan University of Technology),Wuhan 430063,China;2 Shanghai Maritime Safety Administration,Shanghai 200086,China;3 Key Laboratory of High Performance Ship Technology(Wuhan University of Technology),Ministry of Education,Wuhan 430063,China)

In this paper,a deterministic model of ship under keel clearance is presented,which contains a set of random variables caused by ship maneuvering,including the vertical motion induced by wave,the heel during ship turning,the squat in sallow water and hydrographic data of seabed.They are represented by certain distributions,and the ship basic response to typical wave is predicted by a CFD method.Via Monte Carlo simulation,the stochastic uncertainty of a ship under keel clearance is observed,and the grounding probability can support channel construction and vessel traffic safety analysis in restricted waterways.

Monte Carlo;dynamic draft;ship grounding probability

U657.2 U676.1

A

10.3969/j.issn.1007-7294.2014.05.007

1007-7294（2014）05-0532-08

2013-11-23

国家自然科学基金委员会（NSFC）与荷兰科学研究组织（NWO）合作研究项目（51061130548）

任亚磊（1988-），男，武汉理工大学硕士研究生，E-mail：renyalei1988@163.com；

牟军敏（1974-），男，博士，武汉理工大学教授。