基于BPNN的球艏降阻优化模型构建研究

2021-03-22张维英，周俊秋，于博文，于洋

大连理工大学学报 2021年2期

张维英，周俊秋，于博文，于洋

(大连海洋大学航海与船舶工程学院，辽宁大连 116023)

0 引言

船舶快速性是船舶设计中要考虑的主要问题之一，减少船舶航行过程中的阻力，是提高船舶快速性的主要途径.近年来，基于CFD的船舶阻力预报和船型快速性优化取得了显著的成果.叶萌等利用若干设计变量样本点的CFD计算结果，通过建立BP神经网络(BPNN)近似模型成功对载重46 000 t油船进行了型线优化[1].董震鹏使用CFD对某拖网渔船成功进行了球艏形状降阻优化.该文献指出由于对拖网渔船模型采用的是全船模拟法，导致软件计算时间过长，对计算机硬件要求较高，此方法不利于在中小型工作站的普遍推广[2].

对于船型优化中产生的大量船型样本，可根据分阶段采用变精度模型(variable fidelity model，VFM)概念，先后使用低与高两种精度模型对船体进行阻力预报，经过高精度(high fidelity，HF)模型对低精度(low fidelity，LF)模型计算结果的补偿后，得到逼近实际值的计算结果.采用VFM概念进行阻力预报相比于仅用高精度模型进行仿真计算而言，将大幅减少计算量，节约优化成本.Robinson等利用VFM在机翼设计和蝙蝠飞行问题的研究上大幅度降低了耗时[3].Toal也验证了使用VFM比只使用高精度模型更能节约计算成本[4].

求解N-S方程的CFD软件具有良好的精度和可信度，可作为高精度模型的计算工具.高精度模型的计算值在优化过程中不必在每次迭代都进行更新，但低精度模型在每次迭代时都需要进行更新，故寻求一种使用方便且结果精确的低精度模型是分阶段采用VFM概念进行球艏优化的关键.林炜杰使用Holtrop法作为低精度模型对船体进行阻力计算，并成功将VFM应用于船型初步设计[5]，该文献的Holtrop法阻力预报是基于高质量的三维模型进行计算，每一次优化迭代使用Holtrop法进行阻力预报前均需要依据设计点参数值对船模的尺寸进行修改，导致计算效率低、计算成本高，不利于VFM应用于工程实践.本文采用一种将近似技术与Holtrop法相融合求解低精度模型的思路，得到一种利用数学近似函数表达式来快速预报大量设计点的低精度模型.近似技术的出现提供了一种新的优化问题求解方法，其基本思想是依据有限数量的试验数据，构建一个既满足精度要求，计算成本又相对小的数学近似函数表达式用于模拟设计优化问题的输入输出关系，并用于代替实际仿真程序[6].张舒应用BPNN原理，结合试验资料与图谱构造了神经网络阻力预报模型，证明了人工神经网络在阻力预报方面的可行性[7].Volpi等使用随机径向基函数(radial basis functions，RBF)建立动态径向基近似模型，对一艘双体船的阻力和耐波性进行优化，该方法展现出很高的优化精度和效率[8].

本文以一艘参数化试验油船模型球艏的长度与高度快速性优化为例，通过基于拉丁超立方设计法选取的样本点及Holtrop法阻力预报响应值构建基于BPNN的球艏降阻优化低精度模型，并使用该模型对油船进行大量设计点的快速阻力预报，通过分析设计变量对总阻力的影响特性，得到逼近最优解的设计点.

1 标准模型阻力预报研究

本文将一艘75 000 DWT散货船作为标准模型，采用Rhino软件对其进行三维建模，并采用Grasshopper(简称GH)软件的船舶设计插件Nemo 1.1 Beta 的Holtrop法阻力预报功能对75 000 DWT散货船标准模型进行一系列航速下的阻力预报.75 000 DWT散货船主要参数见表1，船模的缩尺比为1∶1.图1展示了该75 000 DWT散货船的三维模型.

图1 75 000 DWT散货船三维模型

表1 75 000 DWT散货船主要参数

Holtrop依据334个模型的试验数据，加上系列64的高速船数据进行回归，使得Holtrop法运用范围扩展至弗劳德数为0.55以上.艾亚法适用范围较广，常应用于中、低速商船，但其统计资料仅代表20世纪40年代以前的船型，故用于近代船型与大型丰满船型时通常误差较大.兰泼凯勒法将船的剩余阻力和摩擦阻力分开计算，剩余阻力由曲线图谱查得，适用于单螺旋桨大型船舶.图2为75 000 DWT散货船的Holtrop法阻力预报值与艾亚法阻力预报值、兰泼凯勒法阻力预报值及试验值的对比[9].

图2 多方法的阻力值比较

比较图2中各方法的阻力预报值，参照各方法的性质，可知Holtrop法的误差与其他两种传统经验公式法相差不大，且船型适用范围更广，Holtrop法可作为低精度模型的求解方法.

2 基于Rhino与GH的船体建模

通过Rhino及其可视化编程插件GH建立本次试验的优化对象——一艘油船三维模型，其基准设计模型的主尺度参数如表2所示.

表2 试验油船模型主尺度参数

试验油船基准设计模型的建立是依据表3所展示的油船船型参数取值范围[10].试验油船基准设计模型如图3所示.

表3 油船船型参数取值范围

图3 油船三维模型

3 神经网络的选用及其应用机理

3.1 神经网络的选用

本文在低精度模型中引入近似技术，提出一种将BPNN作为球艏降阻优化低精度模型的方法.

工程问题中较为常用的近似模型形式有响应面模型(RSM)、Kriging模型和神经网络模型.文献[11]对以上3种近似模型形式进行了研究与对比.由于神经网络模型适用范围广，预测空间相对平滑，拟合效率较高、耗时短，且阻力计算采用神经网络模型精度最高，其性能指标的平均绝对百分比误差(MAPE)较其余近似模型低，故本文选用在船型优化设计领域应用较为成熟的BPNN作为本次优化问题的近似模型形式.

3.2 神经网络的应用机理

只具有一个神经元的单层感知机模型如图4所示，单层感知机功能十分有限.其原理是通过均匀试验设计法选取多组样本点，通过训练神经网络，计算得出满足误差要求的权重，则可以生成数学近似函数表达式用于模拟问题的输入输出关系.

图4 单层感知机模型

设输入为X=(x1x2…xn)T∈Rn，相应权值为W=(w1w2…wn)T∈Rn，则单层感知机的输出y可表示为[12]

(1)

当采集的是高度非线性的样本点时，则需要多层感知机BPNN模型来训练样本点，根据Kolmogorov 定理，只具有一个隐含层的三层BPNN就可以实现对任意精度的逼近.BPNN由输入层、隐含层和输出层组成，通过将输出层的校正误差沿原信号传递通路反传，并逐层调整输出层到隐含层、隐含层到输入层的权值.其模型图如图5所示.

图5 三层BPNN模型图

4 基于BPNN的船舶阻力预报方法

4.1 试验设计方法的选取及其应用

BPNN的学习训练需要依据样本点数据，而均匀的样本点需要利用试验设计方法选取.本文采用具有高灵活性和广泛适用性的拉丁超立方设计法选取构建BPNN的样本点[13].

拉丁超立方抽样是一种从多元参数分布中近似随机抽样的方法，属于分层抽样法.该方法根据样本点数量k的要求，将n维设计空间中每一维设计变量平均划分为k个子空间，所有设计变量的子空间组合形成kn个子空间.再根据每一个设计变量的每一个子空间只出现一次的原则，随机选出k个子空间，从而保证取点的均匀性，最后在选出的每一个子空间中随机取一个样本点，形成k个样本点.

4.2 Holtrop法阻力预报原理

本文选用GH的船舶设计插件Nemo 1.1 Beta软件计算样本点响应值.Nemo软件的阻力预报模块是根据Holtrop在1978～1984年发表的资料为计算依据的.Holtrop法总阻力表示为[14-16]

Rt=Rf(1+k1)+Rapp+Rw+Rb+Rtr+Ra

(2)

式中：Rt是总阻力；Rf是摩擦阻力；1+k1是船型黏性阻力因子；Rapp是附体阻力；Rw是兴波阻力；Rb是球艏附加压阻力；Rtr是艉封板浸水产生的附加压阻力；Ra是模型与实船相关修正因子.

Holtrop已在文献[14-16]中对上式中各阻力的回归计算公式进行了详细说明.因为船型设计初期只对裸船体进行研究，故附体阻力为0.通过对其余各阻力回归计算公式的参数归纳，得到总阻力的一般函数表达式：

Rt=f(v，Shull，Lwl，B，T，，Cp，Lcb，Cm，

Cb，Cw，Abt，Atr，Cstern，Hb，Tf，ie)

(3)

式中：v是航速(kn)；Shull是船体湿表面积；Lwl是设计水线长；B是设计水线宽；T是平均吃水；是型排水体积；Cp是棱形系数；Lcb是浮心纵向位置；Cm是舯横剖面系数；Cb是方形系数；Cw是水线面系数；Abt是艏垂线处球艏横剖面面积；Atr是零航速下的艉封板浸水面积；Cstern是船艉形状系数；Hb是Abt的形心到基线的高度；Tf是艏吃水；ie是半进流角.

孙源等使用基于Holtrop法的Maxsurf软件分析DTMB5415标准模型得出不同航速下球艏纵向位置的合理变化可使总阻力降低1.26%～10.31%[17]，说明改变球艏线性尺寸的方法在降阻方式中有着重要地位.

式(3)中受球艏前端点距艏垂线距离lb、球艏前端点距设计水线距离hb影响的变量有以下4个：船体湿表面积Shull、型排水体积、艏垂线处球艏横剖面面积Abt及其形心到基线的高度Hb，故通过改变设计变量lb、hb的值可改变上述4个变量的值，进而改变由Holtrop法计算的总阻力Rt的值.在优化球艏过程中保持油船总长、型深及型宽等总尺度不变，艏垂线位置随球艏线性尺寸变化而做较小的改变.本文的设计变量lb、hb示意见图6.

图6 球艏几何参数

5 基于BPNN的球艏降阻优化

5.1 球艏降阻优化BPNN的构建流程

本文所构建的球艏降阻优化低精度模型是一种基于一定数量样本点及其响应值所训练的BPNN.图7给出了构建球艏降阻优化BPNN的流程.通过将均匀试验设计法选取的样本点及其响应值送入BPNN进行训练，生成可代替船舶性能计算软件进行快速阻力预报的低精度模型.再通过大量设计点的BPNN总阻力计算结果，预报总阻力随设计变量的变化趋势，分析最优设计点的大致位置，确定最优设计点的近似解，验证球艏降阻优化BPNN作为低精度模型的可行性.

图7 BPNN构建流程图

5.2 第一批样本点的选取及球艏变形

依据文献[10]，只有当航速大于某一数值时，球艏才会显示其降阻作用，该航速值称为球艏的界限速度，并将该界限速度作为是否选用球艏的衡量标准[10]，该试验油船界限速度vk对应的弗劳德数计算公式如下：

Frb=0.644-0.641Cb

(4)

对3 500～7 000 DWT的油船，其服务航速通常为11.5～13 kn.通过将试验油船的Cb代入式(4)中，得出该油船球艏的临界弗劳德数Frb，进而算出该临界弗劳德数对应的界限速度vk=12.26 kn.经查阅船型资料，本文以油船在v=12.5 kn时的降阻为优化目标.

本文分别以GH中试验油船基准设计模型球艏的初始长度lb=4 693 mm和初始高度hb=5 015 mm为两设计变量的基准点.lb的取值范围取以基准点为中心，以1 000 mm为半径的闭区间；hb的取值范围取以基准点为中心，以500 mm 为半径的闭区间.设计空间如下：

3 693 mm≤lb≤5 693 mm

4 515 mm≤hb≤5 515 mm

(5)

基于拉丁超立方设计法在设计空间内选取20个样本点作为BPNN的第一批学习样本点，样本点的二维分布如图8所示，样本点数值见表4.

表4 基准点与部分第一批样本点

图8 第一批样本点的分布

依据20个样本点，在GH中对试验油船模型进行球艏变形，生成20个不同球艏几何参数的三维模型，样本1～5的球艏模型示意如图9所示.

图9 部分球艏变形示意图

5.3 第一批样本点的总阻力计算及相关分析

将依据基准点及第一批样本点所生成的不同尺寸球艏的试验油船依次拾取到GH的船舶设计插件Nemo 1.1 Beta中进行阻力预报.图10展示了Nemo的船体输入模块，在Float区域中设置设计水线的高度，Brep电池用于拾取Rhino中的船体曲面.

图10 Nemo的船体输入模块

输入船体曲面后，Nemo即对船体进行流体静力学计算，图11展示了Nemo的流体静力学计算模块.船体数据经电池左侧Hull_Brep接口输入，通过电池运算后输出于其右侧的HH_Json接口.

图11 Nemo的流体静力学计算模块

在进行阻力预报前，要先定义流体的参数和航速范围.本文设定流体的温度为15 ℃，类别为海水.航速上下限分别为15 kn和12 kn，梯度为0.5 kn.Nemo的流体参数及航速设置模块如图12所示.

图12 Nemo的流体参数及航速设置模块

图13是Nemo的Holtrop法阻力预报模块.

图13 Nemo的Holtrop法阻力预报模块

图14是Nemo软件阻力预报模块全局布置图.含有不同功能的电池模块通过正确的逻辑连线耦合，使数据在各模块间传递，用于计算响应值.

图14 Nemo的阻力预报模块全局布置图

图15展示了油船在球艏尺寸基准点处航速为12.5 kn时GH显示面板的阻力计算结果等数据.

图15 Nemo的阻力数值显示面板

使用船舶设计软件Nemo在试验油船球艏尺寸基准点及20个样本点处测得的模型总阻力如表4所示.

总阻力在设计空间的分布如图16所示.通过颜色标尺观察图16中的散点，可以发现总阻力朝着球艏高度最小值和球艏长度最小值方向减小.

图16 第一批样本点的总阻力空间分布

为了量化变量间的相关程度，本文对第一批样本点的两个自变量与总阻力进行Pearson相关分析.相关系数可以量化变量间的相关程度，反映变量之间的变化趋势.变量x和y的Pearson简单相关系数r定义如下：

(6)

式中：k为样本容量，xj、yj为进行相关分析的两变量所对应的样本值.

数据正态分布是Pearson相关分析的前提之一.通过对3组变量进行正态性检验，得出3组变量正态性检验显著性均大于0.05，可直接进行Pearson相关分析.hb与Rt的相关系数为0.929，Sig值小于0.001，故hb与Rt存在极显著的强正相关关系.lb与Rt的相关系数为0.052，表明lb与Rt存在弱正相关关系，Sig值为0.827(Sig值大于0.05表示变量间相关性不显著)，说明lb与Rt之间的相关性不显著.但根据Pearson简单相关系数的检验t统计量公式：

(7)

可知对任意值的r，样本容量k足够多就能使统计量足够大，通过显著性检验的可能性就高.故本文将通过对BPNN计算的大量设计点变量值进行相关分析，来说明整个设计空间内变量间的相关关系.

观察图16的散点在ZX平面和YZ平面的投影，可发现hb与Rt的散点较lb与Rt更接近一条直线，说明hb与Rt的线性相关性强于lb与Rt.但hb数值变化时为了保证总长等主尺度值不变，从而使艏垂线位置发生较小的改变，进而影响了lb的数值.为进一步分析总阻力的变化规律，并基于变化规律扩充用于训练BPNN的样本点数量，本文将在第一批样本点的阻力最小值点附近开展探究球艏长度与总阻力关系的单因素试验，对球艏长度与总阻力进行Pearson相关分析.

5.4 球艏长度与总阻力的相关分析

在表4样本点的总阻力最小值处，将球艏高度值(hb=4 620 mm)恒定不变，取[4 019，4 419]为球艏长度的研究范围，采用均匀设计试验法生成10组球艏长度试验样本，在GH里对试验油船模型基于10组样本生成10个船体三维模型，并使用Nemo软件对10个模型进行阻力预报.数据见表5.

表5数据所对应的散点图及拟合情况如图17所示，分析图17得出：球艏长度单因素试验数据基本符合线性相关关系.

图17 总阻力的二维分布

表5 球艏长度样本及响应值

通过对两变量进行正态性检验，得出球艏长度和总阻力正态性检验显著性均大于0.05，可直接使用Pearson简单相关系数量化球艏长度与总阻力间的相关程度，经计算得到其相关分析结果：lb与Rt的相关系数为0.631，Sig值为0.05，故在lb的[4 019，4 419]范围内，lb与Rt之间存在显著的强正相关关系.

5.5 第二批样本点的选取及总阻力计算

为了增强BPNN在寻优区间的预报精度，本文在第一批样本点的基础上增添40个样本点.

依据图16中散点的分布趋势和单因素试验得出的变量相关分析结果，判断出最优设计点位于整个设计空间的球艏长度最小值和球艏高度最小值附近.故采用拉丁超立方设计法在该范围内选取40个样本点.第二批样本点的选取范围如下：

3 693 mm≤lb≤4 293 mm

4 515 mm≤hb≤4 815 mm

(8)

使用Nemo软件对选取的第二批样本点进行阻力预报，部分第二批样本点及总阻力如表6所示.

表6 部分第二批样本点

第二批样本点的总阻力在设计空间的分布见图18.

图18 第二批样本点的总阻力空间分布

5.6 球艏降阻优化BPNN的训练

将基准点及第一批和第二批样本点的集合作为输入变量，将其响应值即总阻力作为输出变量，在迈实软件里对球艏降阻优化BPNN进行训练.

图19展示了该BPNN训练过程的训练曲线，从图中可以看出，其均方差σ很快收敛于某极小值处，说明该神经网络的训练达到预期效果.本次试验将该极小值设置为1×10-10，最大迭代数设置为500 000，其意义是当全局误差达到1×10-10时或达到最大迭代数即500 000时，该神经网络的训练终止.

图19 BPNN的训练曲线

通过合理的参数设置，可以得到精度高、误差小的神经网络，进而应用于工程优化.当使用一个隐含层就可以满足求解精度的时候，则建立只含一个隐含层的神经网络，这样有利于神经网络的稳定性，故本神经网络只建立一个隐含层.通过调节该BPNN的隐含层节点个数，确定了该BPNN的隐含层节点为6个，此时训练效果较优.

5.7 球艏降阻优化BPNN的偏差分析

球艏降阻优化BPNN原值与计算值的偏差曲线见图20.

图20 BPNN的偏差曲线

本文采用相对均方根误差(RMSE)和平均绝对百分比误差(MAPE)作为衡量近似模型精度的评价指标，依据第一批与第二批样本点的Nemo软件原值与BPNN计算值而得出的RMSE和MAPE分别为171.6和0.092.

由偏差曲线可以看出，该BPNN计算值能够比较完好地接近于原值，特别是在10～21号样本点之间，计算值与原值几乎重合，在其余样本点上该BPNN计算值也能较好地追踪原值，并反映出总阻力的变化趋势.BPNN的两个评价指标均较小，说明该BPNN对船体总阻力的预报精度较高.该球艏降阻优化BPNN可以在满足精度要求的前提下接替Nemo软件对船体进行总阻力预报，满足成为下一步基于变精度模型的球艏降阻优化研究中的低精度模型的要求.

6 球艏设计变量的局部优化

6.1 第一批设计点的选取及寻优

使用训练完成的球艏降阻优化BPNN在整个设计空间预报第一批由拉丁超立方设计法选取的500个设计点总阻力.BPNN的总阻力预报值空间分布如图21所示.

图21 第一批设计点的总阻力空间分布

观察图21的散点在ZX平面和YZ平面的投影，可以判断第一批设计点的设计变量与总阻力数据基本符合线性相关关系.

通过对3组变量进行正态性检验，得出3组变量正态性检验显著性均小于0.05，故不可直接采用Pearson相关分析，可使用Spearman相关系数进行研究.Spearman相关系数通常被认为是排列后变量之间的Pearson线性相关系数，其用在数据的分布使得Pearson线性相关系数不能用来描述或是用来描述将导致错误的结论时，评价两个统计变量的相关性.

经计算得到其相关分析结果：lb与Rt的相关系数为0.214，Sig值小于0.001，故在整个设计空间内lb与Rt存在极显著的弱正相关关系；hb与Rt的相关系数为0.921，Sig值小于0.001，故在整个设计空间内hb与Rt存在极显著的强正相关关系.

6.2 第二批设计点的选取及寻优

依据整个设计空间内lb、hb与Rt的相关分析结论，通过缩小在设计空间的寻优范围，在存在最优解的局部设计空间：

3 693 mm≤lb≤4 500 mm

4 515 mm≤hb≤5 000 mm

(9)

使用拉丁超立方设计法选取第二批500个设计点，并使用球艏降阻优化BPNN进行阻力预报.第二批设计点总阻力预报值空间分布如图22所示.

观察图22的散点在ZX平面和YZ平面的投影，可以判断第二批设计点的设计变量与总阻力数据基本符合线性相关关系.通过对3组变量进行正态性检验，得出3组变量正态性检验显著性均小于0.05，故不可直接采用Pearson相关分析，可使用Spearman相关系数对第二批500个设计点变量进行相关分析.

图22 第二批设计点的总阻力空间分布

经计算得到其相关分析结果：lb与Rt的相关系数为0.318，Sig值小于0.001，故在整个设计空间内lb与Rt存在极显著的弱正相关关系；hb与Rt的相关系数为0.784，Sig值小于0.001，故在整个设计空间内hb与Rt存在极显著的强正相关关系.该结论与整个设计空间的相关分析结论一致.

比较各设计点及样本点的总阻力，得出总阻力最小值点(lb=3 770 mm，hb=4 638 mm)，该点的总阻力为144 257.6 N，该点是局部优化中逼近最优解的设计点.优化前后的球艏几何形状对比如图23所示.

图23 优化前后球艏变形示意图

7 球艏降阻优化结果分析

对优化结果按船舶阻力理论结合Holtrop法原理进行分析.在Holtrop法各阻力成分中，试验油船优化前后变化最大的分别是兴波阻力Rw、摩擦阻力Rf及球艏附加压阻力Rb.优化后该油船的兴波阻力发生了较明显的减小，摩擦阻力及球艏附加压阻力发生了增大，由于Rw的减少量大于其他阻力成分增量的总和，故总阻力Rt减小，达到了提高船舶快速性的目的.由Holtrop法计算的优化前后两种球艏的总阻力及各阻力成分数值对照如表7所示.