舰船在航行过程中，因受到外加作用会产生晃动，在海浪较大情况下，表现更加明显。舰船的晃动会对需要在舰船进行一系列特殊作业(如：舰载机起降作业中指导舰载机安全起降、导弹发射、风浪中航行的控制等^[1])产生巨大影响，甚至危害生命财产安全。若提前对舰船运动姿态进行预测，可以为决策者提供决策依据，减少安全事故的发生，对航海事业具有重要意义。

目前, 国内外对舰船运动姿态极短期预报非常重视并进行了许多研究。近年来，随着混沌系统理论的快速发展，混沌预报方法被引入到船舶运动姿态的预报中来。文献[2]对于舰船运动姿态序列的混沌特性进行了分析，说明舰船运动姿态序列属于低维混沌序列。

已有研究表明，很多低维混沌时间序列可用二阶Volterra自适应滤波器进行较为精确的自适应预报^{[3, 4]}。文献[5]利用最大Lyapunov指数法对船舶运动姿态时间序列样本进行了混沌性判别，针对船舶运动姿态的非线性和不确定性与混沌特性的紧密关系，给出了船舶运动姿态混沌时间序列的二阶Volterra自适应预报模型。文献[6]提出基于Volterra级数的RLS^[7-10]自适应算法，验证其有较好的收敛性并在预测低维混沌系统取得很好的效果。但是，该方法只是对数据进行一次性处理，没有讨论其对于数据流的预测效果如何。本文在Volterra滤波器与RLS自适应算法的基础上，提出一种基于小波变换的RLS的Volterra级数核估计自适应算法，应用于舰船运动姿态数据流的极短期实时预报研究中。

1 理论基础 1.1 RLS算法原理

RLS算法全称是递推最小二乘算法，它是一种最小二乘自适应横向滤波器的递推算法，它是由n-1时刻滤波器抽头权向量的最小二乘轨迹来递推n时刻权向量的最新估计。

设i时刻模型的输入为x(i), x(i-1), …, x(i-n+1)，ε(i)是标准值y(i)与估计值ŷ(i)之差，这里

$ \hat y(t) = {\mathit{\boldsymbol{\varphi }}^{\rm{T}}}(t)\mathit{\boldsymbol{X}}(t) $

(1)

式中：

$ \mathit{\boldsymbol{\varphi }}(t) = {\left[ {{\varphi _1}(t),{\varphi _2}(t), \cdots ,{\varphi _n}(t)} \right]^{\rm{T}}} $

$ \mathit{\boldsymbol{X}}(i) = {[x(i),x(i - 1), \cdots ,x(i - n + 1)]^{\rm{T}}} $

估计误差为：

$ \varepsilon (i) = y(i) - \hat y(i) = y(i) - {\mathit{\boldsymbol{\varphi }}^{\rm{T}}}(t)\mathit{\boldsymbol{X}}(i) $

设i < 0数据为零，则有：

$ \xi (t) = \sum\limits_{i = 1}^t \beta (t,i){\varepsilon ^2}(i) $

最小二乘的性能指标是使ξ(t)最小，式中β(t, i)代表遗忘因子，0 < β(t, i) < 1, i=1, 2, …, t, 遗忘因子使建模过程对非平稳情况下数据统计特性的变动更具有适应性。通常取指数式因子作为遗忘因子，即

$ \beta (t,i) = {\lambda ^{t - i}}(i = 1,2, \cdots ,t) $

式中λ为趋近于1的正常数，在取指数式因子为遗忘因子时有：

$ \xi (t) = \sum\limits_{i = 1}^t {{\lambda ^{t - i}}} {\varepsilon ^2}(i) $

由əξ/əφ=0可得确定最小二乘意义下最优解$\mathit{\boldsymbol{\hat \varphi }}\left( t \right)$的正规方程

$ \mathit{\boldsymbol{N}}(t)\mathit{\boldsymbol{\hat \varphi }}(t) = \mathit{\boldsymbol{M}}(t) $

(2)

式中：n×n相关阵N(t)和n×1互相关M(t)分别为：

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{N}}(t) = \sum\limits_{i = 1}^t {{\lambda ^{t - i}}} \mathit{\boldsymbol{X}}(i){\mathit{\boldsymbol{X}}^{\rm{T}}}(i)}\\ {\mathit{\boldsymbol{M}}(t) = \sum\limits_{i = 1}^t {{\lambda ^{t - i}}} \mathit{\boldsymbol{X}}(i)\mathit{\boldsymbol{y}}(i)} \end{array} $

(3)

将式(3)中i=t的项单独分开，则有：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{N}}(t) = \lambda \left[ {\sum\limits_{i = 1}^{t - 1} {{\lambda ^{t - i}}} \mathit{\boldsymbol{X}}(i){\mathit{\boldsymbol{X}}^{\rm{T}}}(i)} \right] + \mathit{\boldsymbol{X}}(t){\mathit{\boldsymbol{X}}^{\rm{T}}}(t) = }\\ {\lambda \mathit{\boldsymbol{N}}(t - 1) + \mathit{\boldsymbol{X}}(t){\mathit{\boldsymbol{X}}^{\rm{T}}}(t)} \end{array} $

同理有：

$ \mathit{\boldsymbol{M}}(t) = \lambda \mathit{\boldsymbol{M}}(t - 1) + \mathit{\boldsymbol{X}}(t){\mathit{\boldsymbol{X}}^{\rm{T}}}(t) $

根据(2)可知，计算$\hat \varphi \left( t \right)$需要知道相关矩阵N(t)的逆，利用引理1可使计算矩阵的逆的计算简化(特别是n很大时)，而且在t=1, 2, …过程中是递推式进行。

引理  若A、B为两个n×n维正定阵，且满足：

$ \mathit{\boldsymbol{A}} = {\mathit{\boldsymbol{B}}^{ - 1}} + \mathit{\boldsymbol{C}}{\mathit{\boldsymbol{D}}^{ - 1}}{\mathit{\boldsymbol{C}}^{\rm{T}}} $

式中：D为另一个m×m正定阵; C为n×m阵，则有：

$ {\mathit{\boldsymbol{A}}^{ - 1}} = \mathit{\boldsymbol{B}} - \mathit{\boldsymbol{BC}}{\left( {\mathit{\boldsymbol{D}} + {\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{BC}}} \right)^{ - 1}}{\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{B}} $

令N(t)=A, λN(t-1)= B^-1, X(t)=C, 1=D，经推导可得：

$ \mathit{\boldsymbol{\hat \varphi }}(t) = \mathit{\boldsymbol{\hat \varphi }}(t - 1) + \mathit{\boldsymbol{K}}(t)\mathit{\boldsymbol{a}}(t) $

其中a(t)称作新息，定义为

$ \mathit{\boldsymbol{a}}(t) = \mathit{\boldsymbol{y}}(t) - {\mathit{\boldsymbol{X}}^{\rm{T}}}(t)\mathit{\boldsymbol{\hat \varphi }}(t - 1) $

式中：X^T(t)$\mathit{\boldsymbol{\hat \varphi }}\left( {t - 1} \right)\;$表示根据t-1时刻的${\mathit{\boldsymbol{\hat \varphi }}}$对y(t)所做的“旧的”最小二乘估计，所以a(t)又可以称为先验估计误差。

RLS算法如下：

设定初始值为

$ \mathit{\boldsymbol{P}}(0) = {\delta ^{ - 1}}l,\hat \varphi (t) = 0 $

计算t=1, 2, …时

$ \begin{array}{*{20}{c}} \mathit{\boldsymbol{Z}}(t) = {\lambda ^{ - 1}}\mathit{\boldsymbol{P}}(t - 1)\mathit{\boldsymbol{X}}(t)\\ \mathit{\boldsymbol{K}}(t) = {\left[ {1 + {\mathit{\boldsymbol{X}}^{\rm{T}}}(t)\mathit{\boldsymbol{Z}}(t)} \right]^{ - 1}}Z(t)\\ \mathit{\boldsymbol{a}}(t) = \mathit{\boldsymbol{y}}(t) - {{\mathit{\boldsymbol{\hat \varphi }}}^{\rm{T}}}(t - 1)\mathit{\boldsymbol{X}}(t)\\ \mathit{\boldsymbol{\hat \varphi }}(t) = \mathit{\boldsymbol{\hat \varphi }}(t - 1) + \mathit{\boldsymbol{K}}(t)\mathit{\boldsymbol{a}}(t)\\ \mathit{\boldsymbol{P}}(t) = {\lambda ^{ - 1}}\mathit{\boldsymbol{P}}(t - 1) - \mathit{\boldsymbol{K}}(t){\mathit{\boldsymbol{Z}}^{\rm{T}}}(t) \end{array} $

(4)

定义F_n(t)=λF_n(t-1)+η_n(t)f_n(t)，后验估计误差定义为ε(t)=y(t)-φ^T(t)X(t)，则最小二乘估计的加权误差平方和为：

$ \begin{array}{*{20}{c}} {{\xi _{\min }}(t) = {\xi _y}(t) - {\mathit{\boldsymbol{M}}^{\rm{T}}}(t)\mathit{\boldsymbol{\hat \varphi }}(t)}\\ {{\xi _y}(t) = \sum\limits_{i = 1}^t {{\lambda ^{t - i}}} {y^2}(i) = \lambda {\xi _y}(t - 1) + \mathit{\boldsymbol{y}}(t)} \end{array} $

(5)

由式(4)以及式(5)可得:

$ \begin{array}{*{20}{c}} {{\xi _{\min }}(t) = \lambda {\xi _{\min }}(t - 1) + \mathit{\boldsymbol{a}}(t)\left[ {\mathit{\boldsymbol{y}}(t) - {{\mathit{\boldsymbol{\hat \varphi }}}^T}(t)\mathit{\boldsymbol{X}}(t)} \right] = }\\ {\lambda {\xi _{\min }}(t - 1) + \mathit{\boldsymbol{a}}(t)\varepsilon (t)} \end{array} $

1.2 RLS算法应用于Volterra核估计原理

根据RLS算法的收敛性可知，当输入向量相关矩阵的特征值扩展很大，RLS算法使得收敛速度也得以保证，并且RLS算法通过选取合适的自适应滤波器的系数，来保证输出信号y(t)与期望信号尽可能地匹配。

RLS算法中确定目标函数为

$ {\xi ^d}(t) = \sum\limits_{i = 0}^t {{\lambda ^{t - i}}} {e^2}(i) = \sum\limits_{i = 0}^t {{\lambda ^{t - i}}} {\left[ {\mathit{\boldsymbol{d}}(t) - {\mathit{\boldsymbol{U}}^{\rm{T}}}(t)\mathit{\boldsymbol{H}}(t)} \right]^2} $

(6)

式中：e(i)表示i时刻的输出误差，d(t)表示输入信号为U(i)向量时的期望输出；并且H(t)、U(t)分别表示输入向量及自适应滤波器的系数向量；λ为指数加权因子，应该在0 < λ < 1范围内进行选择。

定义系统的输入矩阵为

$ \mathit{\boldsymbol{P}}(t) = {[\mathit{\boldsymbol{U}}(m),\mathit{\boldsymbol{U}}(m + 1), \cdots ,\mathit{\boldsymbol{U}}(t - 1)]^{\rm{T}}} $

式中，U(t)与式(6)中的U(t)相同，m代表二阶Volterra滤波器的输入项数。则输出为

$ \mathit{\boldsymbol{Y}}(t) = {[y(m),y(m + 1), \cdots ,y(t - 1)]^{\rm{T}}} $

式中y(i)(i=m, m+1, …, t-1)代表系统的实际输出，且满足y(i)=H^T(i)U(i)。

利用传统最小二乘法求解Volterra级数核可得

$ \mathit{\boldsymbol{\hat H}}(t) = {\left[ {{\mathit{\boldsymbol{P}}^{\rm{T}}}(t)\mathit{\boldsymbol{P}}(t)} \right]^{ - 1}}{\mathit{\boldsymbol{P}}^{\rm{T}}}(t)\mathit{\boldsymbol{Y}}(t) $

直接利用Volterra滤波器对舰船运动姿态进行实时在线预报时，随着新数据的不断获取，P(t)、Y(t)维数将不断增大，势必会耗费大量存储空间，为此，使用递推最小二乘法对Volterra级数核进行估计。这样在预报过程中P(t)维数是确定的，可以减少了数据对存储空间的占用，其具体步骤总结如下：

设t+1时刻系统的输入矩阵为P(t+1)，对应输出向量为Y(t+1)，输入输出满足：

$ \mathit{\boldsymbol{P}}(t + 1) = {[\mathit{\boldsymbol{U}}(m),\mathit{\boldsymbol{U}}(m + 1), \cdots ,\mathit{\boldsymbol{U}}(t + 1)]^{\rm{T}}} $

$ \mathit{\boldsymbol{Y}}(t + 1) = {[y(m),y(m + 1), \cdots ,y(t + 1)]^{\rm{T}}} $

记：

$ {\mathit{\boldsymbol{\varphi }}_t} = {\left[ {{\mathit{\boldsymbol{P}}^{\rm{T}}}(t)\mathit{\boldsymbol{P}}(t)} \right]^{ - 1}} $

若增加输入数据x(t+1)后可得：

$ \begin{array}{*{20}{r}} {{\mathit{\boldsymbol{\varphi }}_{t + 1}} = {{\left[ {{\mathit{\boldsymbol{P}}^{\rm{T}}}(t + 1)\mathit{\boldsymbol{P}}(t + 1)} \right]}^{ - 1}} = }\\ {{{\left[ {{\mathit{\boldsymbol{P}}^{ - 1}}(t + 1) + \mathit{\boldsymbol{U}}(t + 1){\mathit{\boldsymbol{U}}^{\rm{T}}}(t + 1)} \right]}^{ - 1}}} \end{array} $

令

$ \mathit{\boldsymbol{\varphi }}_{t + 1}^{ - 1} = \mathit{\boldsymbol{A}},\mathit{\boldsymbol{\varphi }}_t^{ - 1} = \mathit{\boldsymbol{B}},\mathit{\boldsymbol{U}}(t + 1) = C,1 = D $

可得：

$ {\mathit{\boldsymbol{\varphi }}_{t + 1}} = {\left[ {{\mathit{\boldsymbol{P}}^{ - 1}}(t + 1) + \mathit{\boldsymbol{U}}(t + 1){\mathit{\boldsymbol{U}}^{\rm{T}}}(t + 1)} \right]^{ - 1}} $

综上，可对Volterra RLS滤波器的核估计算法总结如下：

依次令t=m+1, m+2, …, m+N, N为样本数据个数。则有：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat H}}(t + 1) = }\\ {\mathit{\boldsymbol{\hat H}}(t) + \mathit{\boldsymbol{K}}(t + 1)\left[ {x(t + 1) - {\mathit{\boldsymbol{U}}^{\rm{T}}}(t + 1)\mathit{\boldsymbol{\hat H}}(t)} \right]} \end{array} $

$ \mathit{\boldsymbol{K}}(t + 1) = {\mathit{\boldsymbol{\varphi }}_t}\mathit{\boldsymbol{U}}(t + 1)/\left[ {1 + {\mathit{\boldsymbol{U}}^{\rm{T}}}(t + 1){\mathit{\boldsymbol{\varphi }}_t}\mathit{\boldsymbol{U}}(t + 1)} \right] $

$ {\mathit{\boldsymbol{\varphi }}_{t + 1}} = {\mathit{\boldsymbol{\varphi }}_t} - \mathit{\boldsymbol{K}}(t + 1){\mathit{\boldsymbol{U}}^{\rm{T}}}(t + 1){\mathit{\boldsymbol{\varphi }}_t} $

式中:K(t+1)表示时变增益矩阵; φ_t与Ĥ(t)的初始值为：

$ {\mathit{\boldsymbol{\varphi }}_m} = {\left[ {{\mathit{\boldsymbol{P}}^{\rm{T}}}(m)\mathit{\boldsymbol{P}}(m)} \right]^{ - 1}} $

$ \mathit{\boldsymbol{\hat H}}(m) = {\left[ {{\mathit{\boldsymbol{P}}^{\rm{T}}}(m)\mathit{\boldsymbol{P}}(m)} \right]^{ - 1}}{\mathit{\boldsymbol{P}}^{\rm{T}}}(m)\mathit{\boldsymbol{Y}}(m) $

1.3 小波去噪细节选择

文献[11]中对小波的定义及小波降噪过程已进行较详细的阐述，本文不再赘述。小波降噪中母小波的选择较为重要。而母小波的选择需要考虑到母小波的数学特性及实际应用情景，以保证精度需求。

其中，母小波的数学特性主要包括正交性、正则性、对称性、高阶消失矩、紧支性^[12]。考虑姿态预测需要母小波具备上述数学特性，故最终选择可以很好满足这些数学特性的Daubechies小波函数，即dbN系列小波，经试验验证，该系列中db6小波在分解层数为6时降噪效果最优。

1.4 基于小波去噪的Volterra-RLS实时预报方法

对于Volterra +RLS自适应预报模型，其预报函数是在拟合混沌序列的基础上获得的，其中，运动姿态序列本身就不可避免的噪声，故经过一系列变换之后，使其得到的预报值包含更多的噪声，使得预报精度不仅要受到模型的影响，还受到噪声影响。同时，Volterra +RLS算法只能对”静”的姿态序列进行单次预报，还不具备对运动姿态数据流进行实时预报的能力。

运动姿态时间序列属于非平稳时间序列，对于分平稳时间序列的降噪方法有很多，但是较为常用的经验模态分解(empirical mode decomposition，EMD)、EMD的优化方法集合经验模态分解(ensemble empirical mode decomposition，EEMD)及小波变换。由于EEMD是对EMD的优化形式，故本文这里主要讨论EEMD和小波降噪方法。EEMD是对序列添加白噪声并且进行多次EMD分解，所以其降噪需要时间较多，本文需要进行实时预测，对于算法的效率要求较高，故对于降噪本文最终选取小波阈值降噪方法。

同时考虑Volterra +RLS自适应预报模型存在的问题，本文最终采用的运动姿态时间序列数据流预测的整体框架如图 1所示。传感器采集到的运动姿态数据数据流进入工作站，首先利用滑动窗口技术进行处理，生成概要数据结构，然后概要数据结构经小波分解分解获得各尺度系数，经阈值处理去除噪声，再经Volterra+RLS算法进行预报，得到最终预报值。

对于概要数据结构通过滑动窗口获得，对于窗口长度的选择，综合考虑算法时间复杂度、预报精度及计算机资源占用等问题后，选用窗口长度为400，预报时长为15 s。

2 仿真

对于舰船运动来说，横摇，又称侧滚角，指浸于水中的物体绕最长延伸方向或波浪入射方向的水平轴的旋转振荡运动，即以船舶重心所在的前后轴线(纵轴线)为中心的回转摇晃。船舶在海上最容易发生横摇且摇摆幅度最大，故本文主要研究横摇运动姿态。

本文将某船模总长12.48 m、船宽1.568 m、设计吃水0.404 m、排水量4.672 t、航速6 kn，在浪向45°、浪向90°工况中20 min横摇运动姿态作为试验数据。

每次选用400个数据作为建模样本，对Volterra+RLS、EEMD+Volterra+RLS及小波+Volterra+RLS等3种预报模型进行预报仿真分析。时序序列实时预报中常用均方根误差RMSE对预报结果进行分析，故本文也以均方根误差RMSE对精度进行评估。

1) 预报精度

先以浪向45°横摇为例，利用3种预报模型进行预报起始点为2 182 s，建模数据长度为400，预报时长为15 s的姿态运动预报，预报结果如图 2。从图 2中可以看出在使用Volterra+RLS对运动姿态序列预报前使用EEMD、小波阈值降噪使得预报精度较Volterra+RLS模型预报精度均有所提高，并且小波阈值降噪效果更好。

为了验证在不同工况下，小波阈值+Volterra+RLS的适用性，再选用90°工况下部分横摇试验数据进行预报仿真，对预报误差记录如表 1所示。

从表 1中可以看出，EEMD+Volterra+RLS和小波+Volterra+RLS对于Volterra+RLS来说预报精度都有所提高，但小波+Volterra+RLS的预报精度更高，这样再次对图 2进行验证。

2) 时间复杂度

根据表 1所示，可以分析得出：

a) 对于相同的横摇运动姿态序列，预报时长相同时Volterra+RLS算法多次进行预报耗时基本稳定，而后两种预报方法耗时会产生一定范围的波动。

b) 由于EEMD需要进行多次EMD分解，并寻找符合条件的IMF分量，所以利用EEMD进行降噪的Volterra+RLS预报算法耗时最长。

c) 通过小波降噪优化的Volterra+RLS预报算法耗时可以发现，小波分解需要的时间远小于EEMD分解所需要的时间，尽管小波降噪耗时比Volterra+RLS要长，但预报15 s需要时间稳定在0.2 s，还是完全可以满足预报需求的。

3 结论

本文提出一种基于小波去噪的Volterra+RLS实时预报方法应用到数据量大、流速快的运动姿态数据流的极短期实时预报中。首先对模型中涉及到的基本理论进行介绍，考虑到噪声对精度影响及基于RLS的Volterra核估计算法只能对横摇运动姿态进行单次预报，不能很好地处理动态的横摇运动姿态数据流等问题，提出使用对于时序数据流常用的滑动窗口方法，对实时数据流进行概要数据结构的获取，然后在概要数据结构上进行小波阈值处理及结合Volterra+RLS方法进行预测。对于算法验证工作，则采用不同工况下横摇运动姿态试验数据对本文提出的模型与EEMD+Volterra+RLS及Volterra+RLS模型进行仿真，并从预报精度和时间复杂度两个方面进行分析讨论，经讨论得出，本文提出的基于小波去噪的Volterra+RLS实时预报方法在预测精度和耗时方面较其他两种方法都具有明显优势。