0.   引　言

智能机舱作为智能船舶的重要组成部分，需要利用船舶机舱系统的设备状态监测参数，对异常情况进行报警，并对运行状态和健康状况进行智能分析评估，从而为系统设备的使用、操控、检修、管理等提供决策支持^[1]。

目前，船舶动力系统的设备状态评估和监控报警一般采用固定阈值法进行判断，其中固定阈值可以根据出厂试验、实船测试和经验分析等予以确定。常用的阈值确定方法包括置信区间阈值法、均值方差标准阈值法^[2]和阈值走廊方法^[3]等。如果阈值设置过高，会使报警灵敏度降低，进而导致报警系统无法正常工作；如果阈值设置过低，则会使虚警率升高，进而导致报警泛滥。此外，支持向量机预测模型已被广泛应用于状态监控^[4]，而船舶动力装置智能诊断系统也被开发并应用于船舶系统状态评估^[5]。然而，上述工作都没有根据设备运行的实际情况进行阈值的自适应更新，这将导致监控系统的误报警以及对设备运行基准与参数的判断偏差，最终无法完成准确的智能控制及辅助决策。

基于上述问题以及智能船舶的实际需求，本文拟提出一种回归支持向量机（support vector regression，SVR）预测模型与滑动时间窗相结合的自适应阈值模型。首先，采用模拟退火算法优化支持向量机模型中的超参数C和g，并对船舶动力系统设备的常规状态特征参数进行建模预测；然后，通过Johnson分布体系将建模残差数据转化为正态数据，并利用滑动时间窗方法计算自适应阈值，用于智能船舶动力系统设备的监控报警和状态评估。该方法可以根据设备运行过程中的实时状态参数确定自适应阈值，以减小阈值带宽并降低虚警率，从而提高船舶动力系统设备状态评估的精确度。

1.   自适应阈值模型

1.1   状态特征参数

作为船舶机舱的重要组成部分，智能船舶动力系统主要包括主推进柴油机、柴油发电机组、推进轴系、燃油系统、润滑油系统和冷却水系统等。动力系统设备的常规状态特征参数包括压力、温度、流量和转速等，然而，对于振动类数据而言，由于其自身频率过高，而本文模型计算的时效性不足，故此类数据不在本文的考虑范围之内。设备的运行状态可以采用单参数法或多参数法进行描述，例如，主推进柴油机的运行状态可以选择排烟温度进行单参数描述，而冷却系统的运行状态则选择压力、压差和温度等多参数进行描述。由于船舶动力系统的工况复杂多变，受其自身固有的不确定性以及因负载、操作和工作环境变化不确定性的影响，所以智能船舶动力系统设备的常规状态特征参数阈值的变化范围较大，如果仍然采用固定阈值进行监控报警和状态评估，将导致较大的评估偏差。为解决这一问题，本文将选择智能船舶动力系统的常规状态特征参数作为样本集，并分析自适应阈值方法。

1.2   回归支持向量机预测模型

回归支持向量机是在支持向量机（support vector machine，SVM）的基础上发展而来，其考虑了样本受退化和工况等多种因素影响的结果，具有较强的鲁棒性，适于解决回归预测问题。本文将选用ε-SVR模型对智能船舶动力系统的状态参数进行预测^[6-7]。

假设给定的训练样本集${{D}}$和函数集${{F}}$为：

式中：${x_i}$，${y_i}$分别为第$i$个输入样本数据和第$i$个输出样本数据；R为实数集；n为样本维度；$f\left( x \right)$为非线性回归函数；w为权向量，${{{w}}^{\rm{T}} }$为权向量的转置；b为偏置量。

回归问题的求解过程即为找到一个函数$f\left( x \right) \in {{F}}$，以使该函数在训练样本$x\left( {x \in {{\bf{R}}^n}} \right)$上的值$f\left( x \right)$与样本期望值$y$的误差小于或等于给定的偏差值$\varepsilon$。

选择权向量${{w}}$和偏置量$b$，将下列规划问题进行最优化处理：

式中：${\xi _i}$，$\xi _i^*$为松弛变量；C>0，为惩罚系数，表示对超出$\varepsilon$之外的样本的惩罚程度。

设定核函数$K\left( {{x_i},{x_j}} \right)$，利用拉格朗日乘子，即可推导出其对偶问题：

式中：${a_i}$，${a_j}$分别为第i， j个变量对应的拉格朗日乘子；$a_i^*$，$a_j^*$为最优的拉格朗日乘子；x_j为训练过程中的第j个样本变量。

求得${a_i}$和$a_i^*$之后，即可解出权向量${{w}}$和偏置量$b$，进而得到非线性回归函数：

本文将在Matlab软件环境下，添加 LIBSVM 工具箱以建立SVR预测模型，其中SVR预测模型算法的具体步骤如下^[8]：

1） 选择船舶动力系统设备的常规状态特征参数监测数据作为样本集，并按时序划分为训练集和预测集。

2） 数据预处理。将训练集和预测集进行归一化处理，将其映射到 [0，1] 区间，其中归一化函数为

式中：$x'$为归一化之后的样本数据；${x_{\min }}$和${x_{\max }}$分别为数据集的最小值和最大值。

3） 确定核函数。核函数的选择必须满足Mercer条件（任何半正定的函数都可作为核函数），本文将选择高斯径向基核函数(radial basis function，RBF)：

式中，$\sigma$为核宽度参数。

4） 选择并优化SVR中的超参数C（惩罚系数）和g（核函数中$\sigma$的数值）。

5） 采用优化之后的超参数 C 和 g 对训练集进行训练。

6） 采用测试样本集对SVR预测模型进行预测与验证。

1.3   模拟退火优化超参数

由1.2节可知，SVR模型的核心问题是惩罚系数C和核函数参数g的选取。本文将采用模拟退火（simulated annealing，SA）算法对$\varepsilon {\rm{ - SVR}}$的超参数C和g进行优化，其中Metropolis 准则是模拟退火算法收敛于全局最优解的关键所在，具体如下：

式中：p为接受新解的概率；x_new为新状态，S(x_new)为新解的能量；x_old为当前状态，S(x_old)为当前解的能量；T为当前温度。

通过迭代方式更新目标函数的最优解时，如果新解的能量小于当前解，则直接接受新解；否则，不直接否定新解，而是以一定的概率接受比当前解更差的解，从而使算法跳离局部最优解，并最终获得全局最优解^[9]。如图1所示，假设初始解在图中A点处，算法沿迭代方向收敛于局部最优解B点，但程序不会在此终止，而是根据Metropolis 准则以一定的概率接受比B点更差的解（C点），从而跳出局部最优解，并最终收敛于全局最优解D点。

图  1  模拟退火算法的寻优示意图

Figure  1.  Optimization diagram of simulated annealing algorithm

下载: 全尺寸图片 幻灯片

SA算法的具体实现步骤如下：

1） 初始化。选择足够大的初始温度值T₀，令T = T₀，然后任意选取初始解S₁作为算法迭代的起点，并设定每个温度下的最大迭代次数 L 和结束温度值T_end。

2） 针对当前温度 T 和迭代次数 k （其中 k =1, 2,···, L ），随机产生一个新解S₂。

3） 计算增量∆f = f (S₂) − f (S₁)，其中f (S₁)，f (S₂)分别为S₁和S₂的评价函数。

4） 如果∆f < 0，则接受S₂作为新的当前解，令S₁= S₂；如果概率 p = exp(−∆f / T)大于0~1之间的随机数rand (0，1)，则仍采用S₂替代S₁成为新的当前解，即令S₁ = S₂；否则，保留当前解S₁。

5） 退温操作。设置降温速度q，以及第k次和第k+1次迭代的温度值T_k，T_k+1，其中T_k+1= q×T_k。如果T_k+1< T_end，则停止迭代，输出当前解作为全局最优解，并终止算法；否则，返回步骤2），重复步骤2）~5）。

算法流程如图2所示。SA算法与初始值的选择无关，根据Metropolis 准则，在搜索过程中可以按照一定的概率接受恶化解，从而跳出某局部最优解；此外，该算法具有渐近收敛性，且结构相对清晰^[10]。因此，基于 SA 算法对$\varepsilon {\rm{ - SVR}}$预测模型的参数进行寻优，可以有效提高模型的预测精度。

图  2  模拟退火算法流程图

Figure  2.  Flow chart of simulated annealing algorithm

下载: 全尺寸图片 幻灯片

1.4   残差数据正态转化方法

根据SA-SVR模型计算的预测数据和真实数据之间的残差，即可得出智能船舶动力系统设备运行状态的评估误差，具体表现为某种未知分布。因此，需对残差数据进行正态转化，在残差符合正态分布假设的前提下计算自适应阈值^[11]。本文将选用Johnson分布体系对建模残差变量X 进行正态转化，建立S_B，S_L，S_U（下标B，L，U分别表示 X 为有界、对数正态和无界）三族分布，具体如表1所示。表1中：h为标准正态分布；$\gamma$，$\eta$和$\delta$，$\lambda$分别为Johnson曲线的位置参数和标度参数。

表  1  Johnson分布系统

Table  1.  Johnson distribution system

分布族	正态变换	参数条件	X条件
SB	$h{\rm{ = }}\gamma {\rm{ + }}\eta \cdot \ln \left( {\dfrac{{X - \delta }}{{\lambda + \delta - X}}} \right)$	$\begin{array}{c}\eta ,\lambda > 0, - \infty < \gamma < \infty , \\- \infty < \delta < \infty\end{array}$	$\delta < X < \delta + \lambda$
SL	$h{\rm{ = }}\gamma {\rm{ + }}\eta \cdot \ln \left( {X - \delta } \right)$	$\begin{array}{c}\eta > 0,-\infty < \gamma <\infty , \\-\infty <\delta <\infty\end{array}$	$X > \delta$
SU	$h{\rm{ = }}\gamma {\rm{ + }}\eta \cdot {\sinh ^{ - 1}}\left( {\dfrac{{X - \delta }}{\lambda }} \right)$	$\begin{array}{c}\eta ,\lambda > 0, - \infty < \gamma < \infty , \\- \infty < \delta < \infty\end{array}$	$- \infty < X < \infty$

下载: 导出CSV 

| 显示表格

利用表 1 中的函数，即可将建模残差变量 X 转化为标准正态分布 h ，具体步骤如下：

1） 按照 Slifker等^[12] 的方法，选定任意正数 z。

2） 利用标准正态分布的概率密度公式$\varPhi \left( \zeta \right) = \displaystyle\int_{ - \infty }^\zeta {\dfrac{1}{{\sqrt {2{\text{π}}} }}{{\rm{e}}^{ - {u^2}/2}}{\rm{d}}u}$，计算4个对称、等距离的标准正态偏差${\zeta }_{J}=\left\{-3{\textit{z}},-{\textit{z}},\;{\textit{z}},3{\textit{z}}\right\}$的概率密度$\varPhi \left( {{\zeta _J}} \right)$，其中J=1，2，3，4；u为随机变量。

3） 利用关系式${i_J} = s*\varPhi \left( {{\zeta _J}} \right) + 0.5$，计算样本的第 ${i_J}$ 顺序观测值，其中s为样本总数。

4） 对样本数据进行线性拟合，求得拟合函数${g} \left( i \right)$，然后计算样本的分位值${X_{{i_J}}} = {g} \left( {{i_J}} \right)$，分位值设定如下：${X_{ - 3{\textit{z}}}} = {X_{{i_1}}}$，${X_{ - {\textit{z}}}} = {X_{{i_2}}}$，${X_{\textit{z}}} = {X_{{i_3}}}$，${X_{3{\textit{z}}}} = {X_{{i_4}}}$。

5） 令参数$m = {X_{3{\textit{z}}}} - {X_{\textit{z}}}$，$n = {X_{ - {\textit{z}}}} - {X_{ - 3{\textit{z}}}}$，$p = {X_{\textit{z}}} - {X_{ - {\textit{z}}}}$，设定分位值比率$QR = mn/{p^2}$。

6） 利用分位值比率$QR$区分Johnson分布体系中的3个分布族，区分准则为：当$QR < 1$时，选择S_B 分布；当$QR = 1$时，选择S_L 分布；当$QR > 1$时，选择S_U 分布。

7） 选定分布形式后，利用下列公式计算3种转换类型的具体参数。

（1） 对于S_B分布：

（2） 对于S_L分布：

（3） 对于S_U 分布：

8） 根据选定的分布形式，利用表1即可计算正态转化之后的h值。

残差数据经过正态转化之后，即可在Matlab中进行KS-检验 (Kolmogorov-Smirnov test)，如果验证的输出结果为0，即表示转化之后的数据符合标准正态分布。残差数据正态转化前后的正态概率图和频数分布直方图分别如图3和图4所示。

图  3  正态概率图

Figure  3.  Normal probability plot

下载: 全尺寸图片 幻灯片

图  4  频数分布直方图

Figure  4.  Frequency distribution histogram

下载: 全尺寸图片 幻灯片

在图3中，采用蓝色“+”标记绘制每个残差数据点的经验概率，采用红色实线参考线连接数据的第1和第3个四分位数，并采用红色虚线参考线将实线延伸到数据的两端。如果所有的样本点都在红色参考线附近，则假设样本服从正态分布；如果样本不服从正态分布，则“+”将构成一条曲线。由图3可以看出，残差数据在进行正态转化之前，正态概率图中的“+”未沿红色参考线分布，构成了一条曲线，而经过正态转化之后，蓝色“+”标记相对接近红色参考线。由图4可以看出，残差数据经过正态转化之后的频数分布直方图非常接近红色的正态分布曲线，验证了利用本文方法进行正态转化的可行性。

1.5   自适应阈值的实现

将残差数据进行正态转化之后，即可根据表1中的公式进行阈值计算，具体步骤如下所示。

1） 根据拉依达准则，确定标准正态分布h的阈值区间$- v \leqslant h \leqslant v$，其中v为最大阈值。

2） 根据选定的分布形式，利用下列公式计算X的阈值。

（1） 对于S_B 分布：

（2） 对于S_L 分布：

（3） 对于S_U分布：

3） 将$X = {X_{\rm{test}}} - {X_{\rm{SVM}}}$代入步骤2）中，得到智能船舶动力系统设备的运行参数阈值，其中${X_{\rm{test}}}$为试验数据，${X_{\rm{SVM}}}$为预测数据。

（1） 对于S_B分布：

（2） 对于S_L分布：

（3） 对于S_U分布：

基于上述阈值计算流程，引入滑动时间窗来实现阈值的自适应。首先，假设存在一个窗口包含s个建模残差数据，对其中的数据进行正态转化；然后，按照上述方法计算其阈值，并求取平均值作为此时的阈值；最后，保证窗口的数据总数不变，逐帧向后滑动窗口，重复上述计算，依次得到每一个时刻的阈值，从而实现阈值的自适应。由于窗口的大小将直接影响阈值的准确性和对异常数据的敏感性，因此需根据实际工况对窗口大小进行相应的调整^{[3, 13]}。

根据状态监测的历史数据建立预测模型之后，即可模拟队列尾插法，利用实时测量的状态监测值不断更新模型的训练样本和预测样本，并重复当前模型的计算过程。此外，基于RBF的$\varepsilon {\rm{ - SVR}}$预测模型在构建时已经考虑了退化和不同工况条件的影响，所以不需要更新模型即可实现实时自适应阈值的计算。

智能船舶动力系统设备自适应阈值模型的建立流程如图5所示。

图  5  自适应阈值模型的流程图

Figure  5.  Flow chart of dynamic adaptive threshold model

下载: 全尺寸图片 幻灯片

2.   自适应阈值模型验证

2.1   方法验证

本文将选取某船舶动力系统设备的正常运行状态特征参数进行方法验证。在参数选取过程中，如果滑动窗口的尺寸过大，即窗口内的数据量过多，将减弱单组数据对计算结果的影响，进而降低模型的灵敏度，那么阈值结果所反映的参数变化会相应地出现延迟；如果滑动窗口的尺寸过小，将增强随机误差的影响，可能导致模型过于灵敏，进而降低评估结果的准确性。经综合考虑，本文将窗口尺寸设定为50，以兼顾评估结果的延迟性要求和模型的灵敏性要求，自适应阈值的计算结果如图6所示。

图  6  自适应阈值模型的方法验证

Figure  6.  Verification of dynamic adaptive threshold model

下载: 全尺寸图片 幻灯片

图6中的红色实线为智能船舶动力系统设备的传统固定阈值线，蓝色实线为实时状态参数曲线，绿色实线为自适应阈值曲线。由图6可以看出，与传统固定阈值相比，自适应阈值更为契合船舶动力系统设备的实时状态参数，其阈值带宽更为紧凑，可以提高阈值的灵敏度。在监控报警方面，可以改善因固定阈值设置过高而导致的设备出现异常情况却不进行报警的缺陷，有利于优化船舶动力系统设备的维修与管理工作。同时，自适应阈值通过提高阈值精度，可以更准确地反映船舶动力系统设备的实际运行情况，从而为智能船舶动力系统设备的状态评估和辅助决策提供更加合理、有效的参考基准。

2.2   异常数据验证

主推进柴油机的运行状态可以采用单参数，即“排烟温度”进行描述。选取某船舶主推进柴油机的排烟温度异常数据作为验证数据，其中实船数据的采集周期为1 s，对所建立的自适应阈值模型进行异常数据验证，结果如图7所示。

图  7  自适应阈值模型的异常数据验证

Figure  7.  Anomaly data verification of dynamic adaptive threshold model

下载: 全尺寸图片 幻灯片

由图7可以看出，与传统固定阈值相比，自适应阈值在主机排烟温度异常升高时，可以提前92 s 成功检测出异常并报警。由此可见，该方法克服了固定阈值对于动力系统设备对异常工况不敏感的缺点，提高了智能船舶监控报警系统的效率，可为系统状态评估和辅助决策提供更精确的数据。

3.   结　语

本文提出了智能船舶动力系统设备的常规状态特征参数的自适应阈值模型。首先，在阈值自适应优化过程中，采用SA算法优化SVR预测模型，实现了对状态特征参数的预测；然后，对预测数据及实际数据计算残差，并结合Johnson分布系统实现了自适应阈值计算。根据实船运行数据的验证结果，作为参考阈值，自适应阈值与传动固定阈值相比其带宽更窄，在设备异常时具有更好的适应性，可以更早地实现故障预警，从而为船舶管理人员提供更充裕的决策时间。此外，自适应阈值与实时监控数据的相关性更紧密，其变化趋势可以作为智能船舶状态评估及寿命预测的参考基准。