基于实测数据的半潜式平台贝叶斯模型修正方法研究

2024-02-26武文华姚伟岸

船舶力学 2024年2期

李松，武文华，2，姚伟岸

（1.大连理工大学工业装备结构分析优化与CAE软件全国重点实验室，辽宁大连 116024；2.大连理工大学宁波研究院，浙江宁波 315000）

0 引言

半潜式平台由于其甲板面积大、稳性好、造价低等优点[1-4]，在深海石油开采中被广泛使用。目前，水动力学仿真和模型实验是海洋平台设计过程中的重要依据和手段，特别是评估极端海况下平台的运动响应。但由于平台实际结构存在的建造误差、平台人员的流动、设备实时维护、货物的吊装等因素，将导致平台结构设计指标发生变化。依据现场监测数据进行平台设计参数修正可以及时掌握平台设计指标的变化规律，有效地提高平台作业的安全等级。

相比于传统的将不确定因素视为确定性偏差的模型修正方法，将待修正参数和修正结果视为统计概率量是对模型修正的一种全新认识。Beck[5]于1998年首次提出了基于贝叶斯推理方法的模型修正及推断的基本思路，随后对求解方法进行改进，提出了一种基于MH 抽样下的自适应Markov Chain Monte Carlo（MCMC）方法，并采用两自由度模型进行验证。由于贝叶斯模型修正方法考虑了先验信息、测试样本等多种不确定因素，在模型修正中得到了广泛的应用。Lam 等[6]利用结构的模态信息结合马尔科夫链蒙特卡洛方法对一个高层办公楼进行了参数修正；Pepi等[7]提出了一种利用环境激励数据进行贝叶斯模型修正的框架，并在钢索桥上验证了该框架的有效性；万华平等[8]利用贝叶斯方法和高斯过程模型进行了斜拉桥的模型修正；Marse 等[9]采用MCMC 方法在三自由度系统中利用频响函数进行参数的修正；Yuen 等[10]在数据不完整的条件下，利用Gibbs 采样实现高维目标求解。为应对大型结构自由度多、计算效率低的问题，Cheung[11]等进行了优化算法的改进。但总体来讲，基于贝叶斯的模型修正主要利用结构的阵型和模态信息进行开展。对于海洋平台这种大型结构，获取其结构的阵型和模态信息需要不同波浪条件下的响应信息，而且不同工况下，系泊条件及压载水会发生相应变化，获取海洋平台的模态信息比较困难。

在求解后验概率分布的方面，马尔科夫链蒙特卡洛（MCMC）方法需要大量的模拟，在复杂结构中应用贝叶斯推理方法仍然是一个挑战，代理模型技术可以有效地解决上述问题。常用的代理模型包括响应面模型、Kriging模型、径向基函数、神经网络、支持向量回归、多变量插值和回归、多项式混沌展开等[12]。Kurtaran 等[13]经过测试选择了三次响应面作为所需的替代函数，并对其进行了工艺参数和制件尺寸的优化设计；Wang 等[14]将Kriging 代理模型与MCMC 结合用于贝叶斯推断，估计了波散射特性的后验概率密度函数。尽管代理模型技术已经在许多工程问题中得到了成功的应用，但目前其在半潜式平台的参数修正方面的应用仍处于较为初步的阶段。

半潜式平台由于建造误差、人员设备等流动行为直接影响平台的水动力学参数，其主要表征为平台的重心变化、平台的排水量变化和回转半径变化，这些变化会直接影响平台响应幅值算子的变化。因此，本文建立了一种基于贝叶斯原理的半潜式平台参数修正方法，来修正上述指标；基于实测的平台六自由度响应特征，根据设计报告中的荷载施加方式将实测的环境荷载作为输入，通过水动力学仿真得到六自由度响应特征；利用代理模型和MCMC方法实现参数后验分布的求解，并将这种方法成功应用于我国南海某半潜式平台的参数修正。首先构建平台的水动力学仿真模型，根据实测海况信息得到不同平台参数的浮体运动响应，然后通过Kriging 代理模型构建待修正参数与响应特征之间的关系，再结合实测响应特征开展基于贝叶斯方法的平台参数修正，最后根据修正的参数对平台多年一遇的响应设计指标进行重新评估。

1 基于贝叶斯推理的半潜式平台参数修正方法

1.1 平台浮体待修正参数和目标函数确定

贝叶斯模型修正中的关键问题之一是如何选取待修正的参数和目标函数。海洋平台长期服役，平台的物理和几何参数（如重量、重心位置等）会发生变化。同时，海洋平台除了受风、浪、流引起的环境荷载外，还要承受由自身运动偏离平衡位置时所产生的静恢复力FS和因自身运动引起的周围海水流体的反作用力，即流体动力FR。其中，海洋平台的浮体运动与静恢复力之间的关系为FS= -CX，其中C代表平台的回复刚度矩阵，可以写为

式中，ρ为海水的密度，g为重力加速度，Aw为水线面面积，xˉf为水线面漂心坐标，∇是排水量，-- ——GMT和分别为横稳性高和纵稳性高。这里除海水的密度、重力加速度外，其他量均与物体的形状、质量、质量分布有关[15]。

流体动力FR由两部分组成，一部分是与运动加速度成正比的附加惯性力，另一部分是与运动速度成正比的阻尼力。按照线性势流理论，流体动力FR可以表示为

而浮体的运动与受力关系服从牛顿第二定律，即

将上述的静恢复力FS和流体动力FR代入式（3），进行整理可得

式中：H(ω)的平方被称作响应幅值算子（Response Amplitude Operator，简称RAO）；C为刚度矩阵，包括静水刚度和系泊刚度，而静水刚度与重心位置、水线面积、浮心位置和质量有关，系泊刚度与系泊长度和系泊材料有关；M为质量矩阵，与重心位置、惯性质量有关；μ为附加质量系数矩阵，与物体水下的形状和运动方向有关；λ为阻尼系数矩阵，与阻尼系数和结构物形状有关。综合以上因素，根据水动力学软件AQWA 可调参数及给出的平台原设计参数，本文选取了平台的横摇回转半径、纵摇回转半径、艏摇回转半径、平台的重心高度和平台的排水量这五个参数作为待修正的参数。

目标函数主要是基于响应的特征。由于本文主要修正与平台浮体有关的参数，根据静水刚度的表达式（1）可知，只有横摇、纵摇和垂荡这三个波频运动具有静恢复力。同时，由于半潜式平台的水平面内回复刚度是由系泊提供的，所以横荡、纵荡和艏摇这三个低频运动主要受系泊影响。因此，本文选取横摇、纵摇和垂荡这三个波频运动的特征作为目标函数。为同时表征运动的剧烈程度和运动的周期特点，选取了时域的统计信息和频域的信息，其中，方差、四分位差、平均峰值等统计量能够反映波频运动的剧烈程度。但是通过对上述的特征进行相关性分析后发现，上述的统计特征之间相关性较强，最终选取了三个自由度的方差作为时域上的统计特征。由于本文选取的运动为波频运动，因此将功率谱分析中的频率成分进行划分，将波浪频率范围内的功率谱密度作为频域上的特征。综上所述，目标函数为三个波频运动的方差和三个波频运动的波频功率谱密度，分别记为横摇方差、纵摇方差、垂荡方差、横摇波频谱密度、纵摇波频谱密度和垂荡波频谱密度。

1.2 参数修正的贝叶斯理论

贝叶斯模型修正的架构是基于贝叶斯理论：

式中，P(θ)为待修正参数的先验分布，P(θ|Y*)是基于观察到的数据更新的参数的后验分布，P(Y*|θ)是似然函数。贝叶斯模型修正是通过先验分布与似然函数更新参数后验概率分布的过程，而似然函数的构建是该过程的核心。似然函数代表不同参数下实测结果发生的概率，通常假设实测和仿真之间符合线性模型[9]：

式中：Y*∈RNm为实测向量，Nm为观测目标值的个数；Y(θ)∈RNm为水动力仿真输出的向量；θ∈RNθ，Nθ表示待修正参数的个数；ε∈RNm表示测试误差向量，用一个均值为零、协方差矩阵为cov的正态分布来模拟。协方差矩阵cov中的对角元素数值大小反映了实测信息的不确定性，在本文中，这个不确定系数的值与传感器的精度有关，具体表现为

式中，ac表示传感器的测量精度，Ns表示实测数据的组数。

在本文中，待修正的参数为横摇回转半径、纵摇回转半径、艏摇回转半径、平台的重心高度和平台的排水量这五个参数。当得到一组实测数据后，提取三个实测波频运动的方差和三个波频运动的波频功率谱密度，将其引入贝叶斯推断模型中，由式（5）获得在给定实测数据下的五个待修正参数的后验概率密度函数。

当先验分布确定时，如果待修正参数的方差不能确定，为了减少先验分布的影响，可选取均匀分布。本文假设待修正的五个参数的先验分布满足变化范围±20%的均匀分布。

1.3 参数代理模型建立

在构建似然关系时，需要构建待修正参数与目标函数之间的关系，但在实际的关系构建过程中，需要模拟不同待修正参数下的响应特征。采用水动力学仿真的方式来确定待修正参数与目标函数之间的关系比较费时。为了提高计算效率，通过代理模型的方式构建待修正参数与目标函数之间的关系。Kriging 代理模型由于其对非线性函数的良好近似能力和独特的误差估计功能，是目前最具代表性的代理模型方法之一[16]。Kriging作为一种半参数化的插值技术，包含了回归部分和非参数部分，所以Kriging模型由两部分组成：多项式和随机分布。

式中：β是回归系数；f(x)是x的多项式，在设计空间中，提供模拟的全局近似，可以是0阶、一阶或二阶多项式：z(x)为随机分布的误差，提供对模拟局部偏差的近似，具有如下统计特性：

式中：n是已知的设计变量的数量；分别是训练样本点xi和xj的第k个分量；Rk(θk,dk)是带有参数θ的相关函数，表征训练样本点之间的空间相关性。所以，Kriging 代理模型将任意响应值都视为服从正态分布的一个随机变量，使模型不限定于某种特殊的形式，而具有更强的灵活性。常用的Rk(θk,dk)有高斯型、指数型、线型、球型和三次型等。

1.4 贝叶斯模型修正的整体技术路线

首先基于待修正参数的联合概率分布，根据抽取参数样本通过水动力学软件AQWA 计算对应参数下的平台运动响应，然后根据运动响应的特征和待修正的参数构建代理模型，通过代理模型和实测的平台运动响应特征构建似然函数，最后通过MCMC 算法进行后验分布的计算。图1 给出了参数修正的整体技术路线。

图1 半潜式平台贝叶斯模型修正过程Fig.1 Bayesian model updating process of semi-submersible platform

具体的计算流程如下：

（1）确定待修正参数后，利用拉丁超立方抽取待修正参数的样本，通过水动力学软件计算出实测荷载下每组待修正参数的响应时程；

（2）利用代理模型构建待修正参数与仿真响应特征之间的关系；

（3）依据代理模型与实测响应特征构建似然函数，并结合参数的先验分布，最终得到待修正参数的后验分布关系；

（4）给定待修正参数的初值θ和建议分布，随机产生第二个待修正值x1，得到α1；

（5）判断α1与0～1之间随机数μ的大小，决定是否保留当前的参数值；

（6）得到关于θ的链，当θ收敛时，θ对应的分布即为后验分布。

2 南海某半潜式平台参数修正过程及结果分析

2.1 半潜式平台现场监测系统

本文基于南海某半潜式平台的现场实测数据开展了贝叶斯参数修正的研究，其中，监测系统的传感器布置情况如图2所示。

图2 半潜式平台监测系统Fig.2 Prototype monitoring system of semi-submersible platform

现场监测数据主要包括风、浪、流的荷载数据和平台六自由度运动的响应数据。相关的测量指标如表1所示。

表1 监测信息的测量范围和精度Tab.1 Range and accuracy of monitoring information

2.2 现场实测数据选取

本文的现场监测数据取自南海某半潜式平台在台风“洛坦”期间的半小时监测信息。相关的环境数据见表2。

表2 台风“洛坦”实测数据Tab.2 Monitoring data of Typhoon Lotan

平台浮体的六自由度响应实测数据如图3所示。

图3 平台浮体六自由度响应实测数据Fig.3 Time histores of 6 DOF response of semi-submersible platform

2.3 待修正参数与目标函数的关系构建

为构建代理模型，需要采用水动力学仿真的方式进行不同参数下平台运动响应的模拟。半潜式平台承受的环境荷载主要包括风、波浪和海流。在水动力学软件AQWA中，将势流理论与莫里森方程相结合，得到平台所受的波浪力；根据风力系数与流力系数得出风力和流力，相关的计算公式如下：

式中，Cwi和Ccu分别为风力系数和流力系数，ρai和ρwa分别为空气和海水的密度，Awi和Acu分别为平台的迎风面积和迎流面积，νwi和νcu分别为风速和流速，Fwi和Fcu分别为风力和流力，Mwi和Mcu分别为风倾力矩和流倾力矩，Lwi和Lcu分别为风倾力臂和流倾力臂，上述单位均为国际单位制。

为了排除其他参数的影响，本文选取的ρai、ρwa、Awi和Acu等参数均来自设计报告。同时，考虑到流体动力受下浮体的形状影响，本文建立的模型也依托平台的设计报告。最终得到的浮体模型和系泊模型如图4所示。

图4 半潜式平台的水动力学模型Fig.4 Dynamical model of semi-submersible platform

本文通过拉丁超立方抽样方法[17-18]抽取100组待修正参数的样本，并代入水动力学软件AQWA进行平台运动响应的计算，然后提取三个波频运动的特征。由于参数较多，因此选取横摇回转半径为30 m，纵摇回转半径为30 m，艏摇回转半径为33 m 时，重心高度与平台排水量分别与波频运动特征之间的代理模型，如图5所示。

图5 参数与运动特征之间的代理模型Fig.5 Surrogate model between parameters and motion characteristics

由图5可知，平台待修正参数与浮体运动响应之间呈现较强的非线性，且横摇的时域和频域特征主要与横摇的回转半径有关，纵摇的时域和频域特征主要与纵摇的回转半径有关，说明了所建立的代理模型的准确性，同时根据表达式（18）给出了代理模型的R2值。

式中，yi表示一系列的观测值表示代理模型给出的预测值，yˉ表示真实值的均值。

由表3可知，预测的R2值均在0.90以上，说明了构建的代理模型的准确性。

表3 参数与运动特征之间的代理模型的R2Tab.3 R2 of the surrogate model between parameters and motion characteristics

2.4 待修正参数的后验分布

通过贝叶斯进行参数修正，可以得到平台的待修正参数的后验分布直方图（图6），表4 中给出了各参数的95%置信区间。

表4 参数的95%置信区间Tab.4 95%confidence interval for the parameters

图6 待修正参数的后验分布Fig.6 Posterior distribution of the parameters to be modified

由结果可知，参数的后验分布为正态分布，其中平台排水量的后验分布的方差最大，这主要是因为横摇运动的方差、纵摇运动方差和垂荡运动的方差对质量比较敏感。因此，当代理模型有一定误差时，由上述三个特征构建的目标函数之间将存在矛盾，导致结果在真实值附近振荡。

将后验概率最大的参数值代入水动力学仿真模型，重新得到修正后的响应。为了具体说明仿真结果与实测结果之间的区别，定义二者之间的相对误差为

式中，S代表水动力学仿真得到的响应特征，M代表实测得到的响应特征。图7 给出了修正前后仿真响应特征与实测响应特征的误差对比。

图7 修正前后的误差对比图Fig.7 Error comparison before and after correction

可以看出，修正后的模型能够较好地还原时程的时域特征和频域特征，在与实测特征的对比上，误差均小于7%。其中，模型修正后垂荡运动成分占比的误差变大，这可能是代理模型的误差造成的。

3 修正后平台浮体的设计指标评估

海洋平台在设计初期，安全评估所用的海洋环境荷载一般有三种[19]，分别是：

A：百年一遇风荷载与其相关联的波浪荷载、流荷载；

B：百年一遇波浪荷载与其相关联的风荷载、流荷载；

C：百年一遇流荷载与其相关联的波浪荷载、风荷载。

其中，南海流花海域海洋环境荷载参数如表5所示。

表5 南海流花海域三种极端海洋环境荷载参数Tab.5 Three kinds of extreme sea environmental load parameters in Liuhua area of the South China Sea

利用上述三种荷载组合用修正后的水动力学模型对平台的运动性能进行重新评估。计算在三种极端荷载均与船艏向成0°、45°、90°、135°、180°条件下半潜式平台的运动响应。将计算得到的三种环境荷载条件下平台六自由度运动响应的最大值和最小值进行统计，具体的统计值如图8所示。

平台在设计时主要关注平台运动响应的最大值，因此，给出了三类海况下平台六自由度响应的最大值，如表6所示。

表6 三类海况下六自由度的最大值Tab.6 Maximum values of 6-DOF response under three design sea conditions

由图8 和表6 可知，结果满足荷载与响应的对应关系。转动自由度的最大值为A 类海况90°作用下的横摇运动，转动自由度的最小值为B 类海况90°作用下的横摇运动。平动自由度的最大值为A 类海况0°作用下的纵荡运动，平动自由度的最小值为A类海况180°作用下的纵荡运动。

整体来讲，横摇、纵摇和垂荡的最大值与最小值主要出现在B类海况下，横荡、纵荡和艏摇的最大值与最小值主要出现在A 类海况下。而A 类海况代表百年一遇风荷载与其相关联的波浪荷载和海流荷载，B 类海况代表百年一遇波浪荷载与其相关联的风荷载和海流荷载，这也说明横摇、纵摇和垂荡主要受波浪荷载影响，横荡、纵荡和艏摇主要受风荷载影响。

由于平台系泊系统的非对称布置，平台沿船艏方向或者垂直于船艏方向的荷载也能引起其他方向的运动，也可以看出对称角度的两种荷载作用在平台上的响应是不同的。但是根据对不同作用方向的荷载作用下的响应分析可知，该半潜式平台抵御与10 号、11 号锚链布置方向一致的荷载能力更强。

设计规范对半潜式钻井平台在海洋环境下的运动响应有一定的要求，一般希望水平运动控制在工作水深的10%以内，垂荡运动不超过±3 m，摇晃小于3°～5°。本文中的半潜式平台的工作水深为310 m 左右，由此可见：三类海况下的横摇和纵摇的极值略超出设计，纵荡在A类和B类海况下的极值略超出设计，垂荡在三类海况下的极值略超出设计值，艏摇在三类海况下的极值均符合设计。