Ship heading model identification based on full rank decomposition least square method
-
摘要: 为了解决标准遗忘因子最小二乘法在线辨识船舶航向模型参数漂移和发散问题，考虑到船舶在实际航行中存在海洋环境扰动和数据欠激励的情况，提出并验证了一种基于满秩分解的递推最小二乘法。用实船数据进行船舶航向模型参数辨识，将辨识结果与标准遗忘因子最小二乘算法、多新息最小二乘法、最小二乘支持向量机的辨识结果进行对比，验证了满秩分解有效降低了在线辨识过程中扰动导致的参数漂移并成功抑制了参数的发散，提升了遗忘因子最小二乘法的辨识精度，减小了最小二乘法对持续数据激励的依赖。Abstract: In order to solve the problem of parameter drift and divergence in the on-line identification of the ship heading model by the forgetting factor least squares method, considering the marine environment disturbance and data under-excitation in actual navigation, we propose a forgetting factor recursive least square algorithm based on full rank decomposition, which uses the ship navigation data to identify the ship heading model parameters, and compare the identification results with the identification results of the standard forgetting factor least squares algorithm, multi-innovation least square algorithm, least square support vector algorithm. The comparison result shows that the full rank decomposition method can effectively reduce the parameter drift caused by the disturbance in the online identification process and successfully suppress divergence of the parameters, improve the accuracy of the forgetting factor least square algorithm and reduce the dependence of the least square method on continuous data excitation.
-
随着数字化和智能化技术的发展，智能船舶发展呈现出日益加快的趋势[1-3]。我国目前处于《智能船舶发展行动计划2019–2021年》提到的三年目标阶段，正在实现智能船1.0向智能船2.0的过渡。自动靠离泊、远程遥控、自主航行等智能船典型应用场景需要精确的船舶航向自适应控制算法[4-7]。国内已产品化的航向控制器与国外先进产品相比，航向精度、能耗、转向快速性和转向平稳性等性能指标仍有较大差距，主要原因有：1) 航向模型参数在线辨识的实时性、收敛性和辨识精度难以同时保证，在线估计时变参数和扰动的能力不够；2) 海洋高频扰动、测量噪声和测量病态引起的航向角速度信号波动处理困难；3) 舵角控制环应对航速和海况变化的自适应能力不足；4) 操舵执行机构性能受限，与国外接近实现连续控制效果有较大差距。
针对上述第一个问题，目前常用的船舶航向模型参数辨识算法主要分为2类：自校正自适应和模型参考自适应[8]。自校正自适应常用的辨识算法包括最小二乘法[9-14]、支持向量机[15-21]、神经网络[22-24]、极大似然估计[25]、遗传算法[26]等。文献[11]采用最小二乘算法辨识船舶操纵模型的参数，并引入新型的学习率，提升了辨识精度。文献[12]针对船舶响应模型参数辨识问题，在最小二乘算法递推迭代过程中将单新息拓展为多新息，大大改进了辨识算法的准确性。文献[13]为了提升辨识过程的快速性，通过构造辨识误差的评价函数，提出了一种动态遗忘因子最小二乘辨识算法。文献[14]在船舶模型辨识最小二乘算法中引入一个非线性的反正切函数来处理新息，有效地提升了算法的辨识精度。文献[18-20]采用支持向量机对船舶航向模型参数进行辨识研究。支持向量机具有全局最优搜索能力，得到的模型泛化能力强，结合优化算法可以提升辨识精度并简化调参。文献[21]调查了船舶操纵运动的水动力模型，将神经网络应用到船舶操纵模型和三自由度模型辨识中。文献[22]针对大型油船的参数辨识问题，提出了基于人工神经网络的参数辨识方法。文献[23]提出了油轮转向数学模型，通过广义椭球基函数神经网络算法辨识转向模型参数。神经网络辨识模型参数对初值没有依赖性，并且可以保证辨识参数的收敛性。文献[25]针对水面无人艇操纵性预报问题，利用极大似然法辨识无人艇参数，采用辨识的结果预报了无人艇的操纵运动。文献[26]利用船舶操纵模型的航向角数据构造目标函数，然后选择遗传算法进行寻优，提升了辨识效率。模型参考自适应[27]设计方法中，采用构造李雅普诺夫函数，并与航向控制器联合设计，可以保证模型参数辨识李雅普诺夫意义下稳定，同时实现航向渐进稳定。
上述方法对于提高模型参数辨识精度、收敛速度和实时性有很多贡献，但缺乏对辨识算法收敛性分析。工程中通常使用满足数据激励条件的回转、Z型、自由航行实验数据来保证模型辨识算法的收敛性，但在自适应控制过程中，持续激励条件往往无法满足，使得参数辨识瞬时收敛性问题成为自校正自适应控制的难点。对于数据长时间欠激励、受扰条件下如何保证或提升在线辨识算法的收敛性问题还没有被解决。
本文针对数据欠激励受扰系统，设计了一种快速、准确、稳定的线性模型在线辨识通用方法，并将该算法应用于船舶航向模型参数在线辨识过程中。将动态降维的思想引入遗忘因子最小二乘法(forgetting factor least square algorithm, FFLS)，提出了一种基于满秩分解的最小二乘法(full rank decomposition least square algorithm, FRDLS)，使得欠激励数据对降维后的辨识模型重新满足激励条件。所提供的算法能够根据死区阈值动态调整满秩分解矩阵，使最小二乘法选择性更新模型参数，限制了扰动导致的参数辨识误差累计。在样本数据长时间欠激励, 海洋环境干扰工况下，将实船数据用于航向模型参数辨识，对比FFLS和FRDLS算法、多新息最小二乘法[12](multi-innovation least square algorithm, MILS)、最小二乘支持向量机[19](least square support vector machines algorithm, LSSVM)辨识船舶操纵模型参数的辨识结果曲线和航向角速度预报误差曲线，可以发现满秩分解最小二乘法能有效抑制数据欠激励和存在海洋环境扰动时参数辨识过程中的漂移和发散问题，并提高参数辨识的精度，同时具有计算量小、实时性强、简单易实现的优点。
1. 船舶航向模型和基本辨识算法
远洋航行的船舶多工作在航向保持模式下，该工况下，船舶自动舵通过给定小舵角抵抗海洋环境干扰，使实际航向保持在期望航向误差带内，辨识准确的航向模型参数是实现高精度航向控制的前提。
1.1 船舶航向模型
船舶操纵运动模型可以采用一阶线性Nomoto模型[13]：
$$ \left\{ {\begin{split} &{T\dot {\boldsymbol{r}} + {\boldsymbol{r}} = K{\boldsymbol{\delta}} } \\ &{\dot {\boldsymbol{\phi}} = {\boldsymbol{r}}} \end{split}} \right. $$ (1) 式中：
${\boldsymbol{r}}$ 为船体艏摇角速度；$\dot {\boldsymbol{r}}$ 为艏摇角加速度；${\boldsymbol{\phi}}$ 为船舶艏向角；${\boldsymbol{\delta}}$ 为舵角；$ K $ 、$ T $ 为船舶回转能力和回转惯性参数，受到船舶载荷条件和航行速度的影响。最小二乘改进模型辨识算法非常依赖于经典状态空间模型参数，实际工程中考虑到外界扰动和模型摄动对参数辨识的不利影响，此处引入扰动项对Nomoto模型进行改进：$$ \left\{ {\begin{split} &{T\dot {\boldsymbol{r}} + {\boldsymbol{r}} = K({\boldsymbol{\delta}} + {{\boldsymbol{\delta}} _d})} \\ &{\dot {\boldsymbol{\phi}} = {\boldsymbol{r}}} \end{split}} \right. $$ (2) 式(2)中
${{\boldsymbol{\delta}} _d}$ 为干扰力矩等效舵角，包含了舵的零位误差、船舶航向模型非线性、测量噪声中的非零均值部分、海洋环境扰动力矩等。为辨识做准备，将改进Nomoto模型整理为标准最小二乘形式。$$ \dot {\boldsymbol{r}} = \frac{1}{T}{\boldsymbol{r}} + \frac{K}{T}{\boldsymbol{\delta}} + \frac{K}{T}{{\boldsymbol{\delta}} _d} $$ (3) 为了加快参数辨识收敛速度，引入归一化因子
$ {\alpha _{\dot r}} $ 和${\alpha _r}$ ，辨识前通常使用归一化因子使数据向量各维度具有相同幅值范围，令：$$ \left\{ {\begin{split} &{\dot {\boldsymbol{r}}' = {\alpha _{\dot r}}{\alpha _r}\dot {\boldsymbol{r}}} \\ & {{\boldsymbol{r}}' = {\alpha _r}{\boldsymbol{r}}{\text{ }}} \\ &{K' = {\alpha _r}K{\text{ }}} \\ &{T' = T/{\alpha _{\dot r}}} \end{split}} \right. $$ (4) 选取输出量
${{\boldsymbol{z}}_k}$ 、参数向量${{\boldsymbol{\theta}} _k}$ 、数据向量${{\boldsymbol{h}}_k}$ ，下标$k$ 代表第$k$ 时刻：$$ \left\{ \begin{split} &{{\boldsymbol{z}}_k} = \dot {\boldsymbol{r}}'\\ &{{\boldsymbol{\theta}} _k} = {\left( {\frac{1}{{T'}}\frac{{K'}}{{T'}}\frac{{K'}}{{T'}}{{\boldsymbol{\delta}} _d}} \right)^{\rm{T}}}\\ &{{\boldsymbol{h}}_k} = {\left( { - {\boldsymbol{r}}'{\boldsymbol{\delta}} } \right)^{\rm{T}}} \end{split} \right. $$ (5) 整理得最小二乘测量方程为
$$ {{\boldsymbol{z}}_k} = {{\boldsymbol{h}}_k}{{\boldsymbol{\theta}} _k} + {{\boldsymbol{w}}_k} $$ (6) 式中
${{\boldsymbol{w}}_k}$ 为$k$ 时刻测量噪声。1.2 遗忘因子最小二乘法
递推最小二乘算法因其简单、实时性强的优点，广泛应用于模型参数在线辨识。设给定数据集
$\{ {{\boldsymbol{h}}_k},{{\boldsymbol{z}}_k}\} _{k = 1}^N$ ，其中，${{\boldsymbol{h}}_k} = {\left[ {{{\boldsymbol{h}}_{k1}} {{\boldsymbol{h}}_{k2}} \cdots {{\boldsymbol{h}}_{kn}}} \right]^{\rm{T}}}$ 为一次测量数据向量，令数据矩阵为${\boldsymbol{H}} = {[{{\boldsymbol{h}}_1}{{\boldsymbol{h}}_2} \cdots {{\boldsymbol{h}}_N}]^{\rm{T}}} \in {{\bf{R}}^{N \times n}}$ ，输出矩阵为${\boldsymbol{Z}} = {[{{\boldsymbol{z}}_1}{{\boldsymbol{z}}_2} \cdots {{\boldsymbol{z}}_N}]^{\rm{T}}}$ ，则线性回归方程可以表示为$$ {\boldsymbol{Z}} = {\boldsymbol{H\theta}} + {\boldsymbol{w}} $$ (7) 式中：
${\boldsymbol{\theta}} = {[{\theta _1}{\theta _2} \cdots {\theta _n}]^{\rm{T}}}$ 为待估计系数向量；${\boldsymbol{w}} = [ {w_1} {w_2} \cdots $ $ {w_N}]^{\rm{T}}$ 为$N$ 次测量的模型拟合误差向量；${\boldsymbol{\theta}} = [{\theta _1} {\theta _2} \cdots {\theta _n}]^{\rm{T}}$ 为待辨识参数数量。最小二乘法(least square algorithm, LS)采用误差向量二范数为损失函数：$$ {\boldsymbol{J}}(\hat \theta ) = \frac{1}{2}{{\boldsymbol{w}}^{\rm{T}}}{\boldsymbol{w}} = \frac{1}{2}{({\boldsymbol{Z}} - {\boldsymbol{H}}\hat {\boldsymbol{\theta}} )^{\rm{T}}}({\boldsymbol{Z}} - {\boldsymbol{H}}\hat {\boldsymbol{\theta}} ) $$ (8) 令式(8)对
$\hat {\boldsymbol{\theta}}$ 的偏导数等于零，得：$$ {{\boldsymbol{H}}^{\rm{T}}}{\boldsymbol{Z}} = {{\boldsymbol{H}}^{\rm{T}}}{\boldsymbol{H}}\hat {\boldsymbol{\theta}} $$ (9) 矩阵
${\boldsymbol{H}}$ 的秩${\rm Rank}({\boldsymbol{H}}) = n$ 时，${{\boldsymbol{H}}^{\rm{T}}}{\boldsymbol{H}}$ 可逆，由式(9)可得${\boldsymbol{\theta}}$ 的最小二乘估计：$$ \hat {\boldsymbol{\theta}} = {({{\boldsymbol{H}}^{\rm{T}}}{\boldsymbol{H}})^{-1}}{{\boldsymbol{H}}^{\rm{T}}}{\boldsymbol{Z}} $$ (10) 为了将最小二乘转化为递推形式，可令
$k$ 时刻数据矩阵、测量矩阵、方差矩阵和辨识结果矩阵分别为$$ \left\{ \begin{split} &{{\boldsymbol{H}}_k} = {[{{\boldsymbol{h}}_1} {{\boldsymbol{h}}_2} \cdots {{\boldsymbol{h}}_k}]^{\rm{T}}} \\ &{{\boldsymbol{Z}}_k} = {[{{\boldsymbol{z}}_1} {{\boldsymbol{z}}_2} \cdots {{\boldsymbol{z}}_k}]^{\rm{T}}} \\ &{{\boldsymbol{P}}_k} = {({\boldsymbol{H}}_k^{\rm{T}}{{\boldsymbol{H}}_k})^{ - 1}} \\ &{{\hat {\boldsymbol{\theta}} }_k} = {({\boldsymbol{H}}_k^{\rm{T}}{{\boldsymbol{H}}_k})^{ - 1}}{\boldsymbol{H}}_k^{\rm{T}}{{\boldsymbol{Z}}_k} \end{split} \right. $$ (11) 将
${{\boldsymbol{H}}_{k + 1}} = {[{\boldsymbol{H}}_k^{\rm{T}}{{\boldsymbol{h}}_{k + 1}}]^{\rm{T}}}$ 、${{\boldsymbol{Z}}_{k + 1}} = {[{\boldsymbol{Z}}_k^{\rm{T}}{{\boldsymbol{z}}_{k + 1}}]^{\rm{T}}}$ 代入${{\boldsymbol{P}}_{k + 1}}$ 和${\hat {\boldsymbol{\theta}} _{k + 1}}$ 中，令${{\boldsymbol{K}}_{k + 1}} = {{\boldsymbol{P}}_{k + 1}}{h_{k + 1}}$ ，根据矩阵分块运算可得递推最小二乘法：$$ \left\{ \begin{gathered} {{\boldsymbol{P}}_{k + 1}} = {{\boldsymbol{P}}_k} - {{\boldsymbol{P}}_k}{{\boldsymbol{h}}_{k + 1}}{[{\boldsymbol{h}}_{k + 1}^{\rm{T}}{{\boldsymbol{P}}_k}{{\boldsymbol{h}}_{k + 1}} + 1]^{ - 1}}{\boldsymbol{h}}_{k + 1}^{\rm{T}}{{\boldsymbol{P}}_k} \hfill \\ {K_{k + 1}} = {{\boldsymbol{P}}_{k + 1}}{{\boldsymbol{h}}_{k + 1}} \hfill \\ {{\hat {\boldsymbol{\theta}} }_{k + 1}} = {{\hat {\boldsymbol{\theta}} }_k} + {{\boldsymbol{K}}_{k + 1}}[{{\boldsymbol{z}}_{k + 1}} - {\boldsymbol{h}}_{k + 1}^{\rm{T}}{{\hat {\boldsymbol{\theta}} }_k}] \hfill \\ \end{gathered} \right. $$ (12) 为避免随着数据增多，新数据作用逐渐消失，引入遗忘因子
${\beta _1}(0.9 < {\beta _1} < 1)$ 得：$$ \left\{ \begin{gathered} {{\boldsymbol{P}}_{k + 1}} = ({{\boldsymbol{P}}_k} - {{\boldsymbol{P}}_k}{{\boldsymbol{h}}_{k + 1}}{[{\boldsymbol{h}}_{k + 1}^{\rm{T}}{{\boldsymbol{P}}_k}{{\boldsymbol{h}}_{k + 1}} + 1]^{ - 1}}{\boldsymbol{h}}_{k + 1}^{\rm{T}}{{\boldsymbol{P}}_k})/{\beta _1} \hfill \\ {{\boldsymbol{K}}_{k + 1}} = {{\boldsymbol{P}}_{k + 1}}{{\boldsymbol{h}}_{k + 1}} \hfill \\ {{\hat {\boldsymbol{\theta}} }_{k + 1}} = {{\hat {\boldsymbol{\theta}} }_k} + {{\boldsymbol{K}}_{k + 1}}[{{\boldsymbol{z}}_{k + 1}} - {\boldsymbol{h}}_{k + 1}^{\rm{T}}{{\hat {\boldsymbol{\theta}} }_k}] \hfill \\ \end{gathered} \right. $$ (13) 2. 基于满秩分解的最小二乘算法设计
当式(8)中矩阵
${\boldsymbol{H}}$ 的秩$r = {\rm{ Rank}}({\boldsymbol{H}}{)} < n$ 时，${{\boldsymbol{H}}^{\rm{T}}}{\boldsymbol{H}}$ 奇异，不可逆，此时存在正交矩阵${\boldsymbol{V}}$ ，使得：$$ \left\{ \begin{gathered} {\boldsymbol{H}} = [{{\boldsymbol{H}}^*}0]{\boldsymbol{V}} = [{{\boldsymbol{H}}^*}0]\left[ \begin{gathered} {{\boldsymbol{V}}_1} \hfill \\ {{\boldsymbol{V}}_2} \hfill \\ \end{gathered} \right] = {{\boldsymbol{H}}^*}{{\boldsymbol{V}}_1} \hfill \\ {{\boldsymbol{H}}^*} = {\boldsymbol{HV}}_1^{\rm{T}} \hfill \\ \end{gathered} \right. $$ (14) 上标T表示矩阵转置，令
${\hat {\boldsymbol{\theta}} ^*} = {{\boldsymbol{V}}_1}\hat {\boldsymbol{\theta}} \in {{\bf{R}}^r}$ ，整理得：$$ {\hat {\boldsymbol{\theta}} ^*} = {({{\boldsymbol{H}}^{*{\rm{T}}}}{{\boldsymbol{H}}^*})^{ - 1}}{{\boldsymbol{H}}^{*\rm{T}}}{\boldsymbol{Z}} $$ (15) 令
$\hat {\boldsymbol{\theta}} ' = {{\boldsymbol{V}}_2}\hat {\boldsymbol{\theta}} \in {{\bf{R}}^{n - r}}$ ，则有：$$ \left\{ \begin{split} &\left[ \begin{array}{l} {{\hat {\boldsymbol{\theta}} }^*}\\ \hat {\boldsymbol{\theta}} ' \end{array} \right] = \left[ \begin{array}{l} {{\boldsymbol{V}}_1}\\ {{\boldsymbol{V}}_2} \end{array} \right]\hat {\boldsymbol{\theta}} = {\boldsymbol{V}}\hat {\boldsymbol{\theta}} \\ &\hat {\boldsymbol{\theta}} = {{\boldsymbol{V}}^{\rm{H}}}\left[ \begin{array}{l} {{\hat {\boldsymbol{\theta}} }^*}\\ \hat {\boldsymbol{\theta}} ' \end{array} \right] \end{split} \right. $$ (16) 任取
$\hat {\boldsymbol{\theta}} ' \in {{\bf{R}}^{n - r}}$ 和${\hat {\boldsymbol{\theta}} ^*}$ 满足式(15)，可根据式(16)构造满足式(9)的$\hat {\boldsymbol{\theta}}$ ，$\hat {\boldsymbol{\theta}} $ 为${\boldsymbol{\theta}}$ 的最小二乘估计，显然${\boldsymbol{\theta}}$ 的最小二乘估计$\hat{\boldsymbol{\theta}}$ 不唯一，且构成$n - r$ 维线性空间。为了使估计结果唯一，可选取$\hat {\boldsymbol{\theta}} ' = {{\boldsymbol{V}}_2}{{\boldsymbol{\theta}} _0}$ ，${{\boldsymbol{\theta}} _0}$ 为${\boldsymbol{\theta}}$ 的先验初值或上一时刻递推值。主成分${\hat {\boldsymbol{\theta}} ^*}$ 最小二乘解唯一，可以写成递推形式，即$$ \left\{ \begin{gathered} {\boldsymbol{P}}_{k + 1}^* = ({\boldsymbol{P}}_k^* - {\boldsymbol{P}}_k^*{\boldsymbol{h}}_{k + 1}^*{[{\boldsymbol{h}}_{k + 1}^{*{\rm{T}}}{\boldsymbol{P}}_k^*{\boldsymbol{h}}_{k + 1}^* + 1]^{ - 1}}{\boldsymbol{h}}_{k + 1}^{*{\rm{T}}}{\boldsymbol{P}}_k^*)/{\beta _2} \hfill \\ {\boldsymbol{K}}_{k + 1}^* = {\boldsymbol{P}}_{k + 1}^*{\boldsymbol{h}}_{k + 1}^* \hfill \\ \hat {\boldsymbol{\theta}} _{k + 1}^* = \hat {\boldsymbol{\theta}} _k^* + {\boldsymbol{K}}_{k + 1}^*[{{\boldsymbol{z}}_{k + 1}} - {\boldsymbol{h}}_{k + 1}^{*{\rm{T}}}\hat {\boldsymbol{\theta}} _k^*] \hfill \\ \end{gathered} \right. $$ (17) 式中：
$$ \left\{ \begin{gathered} {\boldsymbol{P}}_k^* = {{\boldsymbol{V}}_1}{{\boldsymbol{P}}_k}{\boldsymbol{V}}_1^{\rm{T}} \hfill \\ \hat {\boldsymbol{\theta}} _k^* = {{\boldsymbol{V}}_1}{{\hat {\boldsymbol{\theta}} }_k} \hfill \\ \end{gathered} \right. $$ (18) 正交投影变换矩阵
${\boldsymbol{V}}$ 的选取不唯一，本文结合递推最小二乘法给出${\boldsymbol{V}}$ 的一种求解方法。设第$k$ 次测量数据$\{ {{\boldsymbol{h}}_k},{{\boldsymbol{z}}_k}\}$ 对应的变换阵为${{\boldsymbol{V}}_k} = [{{\boldsymbol{V}}_{k1}}{{\boldsymbol{V}}_{k2}}]。$ $$ \left\{ \begin{gathered} {\boldsymbol{V}}_{k1}^p = [{\boldsymbol{V}}_{k1}^{p - 1}{{\boldsymbol{e}}_i}] , {{\boldsymbol{h}}_{ki}} > {{\boldsymbol{h}}_s} \hfill \\ {\boldsymbol{V}}_{k2}^q = [{\boldsymbol{V}}_{k2}^{q - 1}{{\boldsymbol{e}}_i}] , {{\boldsymbol{h}}_{ki}} \leqslant {{\boldsymbol{h}}_s} \hfill \\ \end{gathered} \right. $$ (19) 式中：
${{\boldsymbol{e}}_i}$ 为第$i$ 行为1，其余行为0的$n$ 维列向量；${{\boldsymbol{h}}_s}$ 为死区阈值，用于削弱干扰的影响；${{\boldsymbol{h}}_{ki}}$ 为第$k$ 次测量数据向量的第$i$ 行数据，当$i = 1,2,\cdots,n$ 时，有$p = 1,2,\cdots,r$ 或$q = 1,2,\cdots,n - r$ ，最终可得：$$ \left\{ \begin{gathered} {{\boldsymbol{V}}_{k1}} = {\boldsymbol{V}}_{k1}^r \hfill \\ {{\boldsymbol{V}}_{k2}} = {\boldsymbol{V}}_{k2}^{n - r} \hfill \\ \end{gathered} \right. $$ (20) 满秩分解最小二乘法算法实现步骤为
1) 初始化式(13)中方差阵
${{\boldsymbol{P}}_k}$ 、参数向量${{\boldsymbol{\theta}}_k}$ ；初始化上一时刻正交变换阵${{\boldsymbol{V}}_{k - 1}}$ 为$n$ 阶单位阵；初始化满秩分解后方差阵${\boldsymbol{P}}_k^*$ 等于${{\boldsymbol{P}}_k}$ ，满秩分解后参数向量${\boldsymbol{\theta}} _k^*$ 等于${{\boldsymbol{\theta}} _k}$ ；2) 通过测量值向量
${{\boldsymbol{h}}_k}$ 根据式(19)、(20)求取正交变换矩阵${{\boldsymbol{V}}_k}$ ；3) 若
${{\boldsymbol{V}}_k}$ 不等于${{\boldsymbol{V}}_{k - 1}}$ ，根据式(18)求${\boldsymbol{P}}_k^*$ 和${\boldsymbol{\theta}} _k^*$ ，否则跳到4)；4) 根据式(13)求
${{\boldsymbol{P}}_{k + 1}}$ ；根据式(17)求${\boldsymbol{P}}_{k + 1}^*$ 和${\boldsymbol{\theta}} _{k + 1}^*$ ；根据式(16)由${\hat{\boldsymbol{\theta}}_k}$ 求$\hat {\boldsymbol{\theta}} {'_k}$ ，令$\hat {\boldsymbol{\theta}} {'_{k + 1}}$ 等于$\hat{\boldsymbol{\theta}} {'_k}$ ；5) 由
${\boldsymbol{\theta}}_{k + 1}^*$ 、$\hat {\boldsymbol{\theta}} {'_{k + 1}}$ 和${{\boldsymbol{V}}_k}$ 根据式(16)求${\hat {\boldsymbol{\theta}} _{k + 1}}$ ；6) 循环2)~5)可实现参数在线辨识。
3. 实验验证
用于验证辨识算法的实验数据为三级海况下，一艘货船的自动舵数据记录仪从大连到新加坡航行过程中记录获得，系统采样频率为50 Hz。艏摇角速度
${\boldsymbol{r}}$ 、舵角${\boldsymbol{\delta}}$ 和航向${\boldsymbol{\psi}}$ 数据如图1所示，角速度信号测量获得后，经过IIR低通滤波器滤除高频测量噪声后用于模型辨识。
图 1 实验数据曲线Fig. 1 Experiment data curves
下载:
全尺寸图片
3.1 船舶航向模型参数辨识
分别采用LS、FFLS、MILS、LSSVM和FRDLS算法，基于图1中的实验数据对式(5)中的模型参数进行辨识，对比辨识结果，分析在辨识数据充分激励和欠激励工况下，辨识算法的辨识精度、收敛性和实时性。
MILS算法新息数选为10，FRDLS算法与FFLS辨识算法遗忘因子均为
$\beta = 0.999\;94$ ，FRDLS算法阈值参数$ {h_s} = {\text{1}} $ ，LSSVM算法选择惩罚因子C为20，数据滑动窗口长度为500 s，K、T的初值选为0.2和250，扰动${\boldsymbol\delta _d}$ 初值为0，式(4)中归一化系数分别为${\alpha _{\dot r}} = 25$ 和${\alpha _r} = 50$ 。对比5种辨识算法，参数辨识结果曲线如图2、3所示。
图 2 K辨识结果曲线Fig. 2 K-coefficient curves下载:
全尺寸图片
图 3 T辨识结果曲线Fig. 3 T-coefficient curves下载:
全尺寸图片
从图2、3可以看出，FFLS、MILS算法辨识K、T参数过程中，在200 s和3 000 s附近均出现参数发散，LSSVM算法数据滑动窗口进入1200秒前后的船舶两次转向数据段时，K、T参数辨识收敛，滑动窗口进入1 800s之后航向保持数据段，辨识参数开始发散；FRDLS算法辨识K、T参数和扰动过程相对平稳，在船舶转向数据段和长时间航向保持数据段均未出现参数发散。辨识算法引入满秩分解后提高了遗忘因子最小二乘法参数辨识过程的鲁棒性。参数辨识结果见表1。
表 1 参数辨识结果Table 1 Parameters identification results算法 $ K $ $T/{\text{s} }$ ${\delta _d}$/(°) LS 0.012 59.5 −0.23 FFLS 0.051 93.2 −0.31 MILS 0.059 116.8 −0.20 LSSVM 0.055 94.3 −0.33 FRDLS 0.038 70.9 −0.36 3.2 船舶运动预报
为了对比上述算法辨识船舶航向模型参数的精度，可以将图2~图4中参数辨识结果代入式(3)中，以舵角数据为输入，对式(3)微分方程进行求解，从而预报船舶艏摇角速度。将预报结果与实验数据作差，可得到预报误差：
$$ {e_i} = {r_i} - {r_{ip}} $$ (21) 式中
${r_i}$ 和${r_{ip}}$ 分别为第$i$ 时刻航向角速度测量值和预报值。预报误差越小，则辨识精度越高。为了保证航向角速度预报精度，此处采用四阶龙格库塔法求解航向模型微分方程。
图 4${\delta _d} $ 辨识结果曲线Fig. 4${\delta _d}$ -coefficient curves下载:
全尺寸图片
$$ \left\{ \begin{gathered} \dot {\boldsymbol{x}} = {\boldsymbol{f}}({\boldsymbol{x}},t) \hfill \\ {{\boldsymbol{K}}_1} = {\boldsymbol{f}}({\boldsymbol{x}},t) \hfill \\ {{\boldsymbol{K}}_2} = {\boldsymbol{f}}({\boldsymbol{x}} + 0.5{T_s}{K_1},t + 0.5{T_s}) \hfill \\ {{\boldsymbol{K}}_3} = {\boldsymbol{f}}({\boldsymbol{x}} + 0.5{T_s}{K_2},t + 0.5{T_s}) \hfill \\ {{\boldsymbol{K}}_4} = {\boldsymbol{f}}({\boldsymbol{x}} + {T_s}{K_3},t + {T_s}) \hfill \\ {{\boldsymbol{ x}}(t + {T_s})} = {\boldsymbol {x}}(t) + {T_s}({K_1} + 2{K_2} + 2{K_3} + {K_4})/6 \hfill \\ \end{gathered} \right. $$ (22) 式中：
${\boldsymbol{x}} = {[\boldsymbol r{\text{ }}{\boldsymbol{\psi}} ]^{\rm{T}}}$ ；${\boldsymbol{f}}({\boldsymbol{x}},t) = {[ (- {\boldsymbol{r}} + K{\boldsymbol{\delta}} {\text{ }})/T {\boldsymbol{r}}]^{\rm{T}}}$ ；${T_s}$ 为系统采样周期，此处取${T_s} = 0.02\;{\rm{s}}$ 。角速度预报曲线和预报误差曲线如图5、6。
图 5 角速度预报曲线Fig. 5 Yaw prediction curve下载:
全尺寸图片
图 6 角速度预报误差曲线Fig. 6 Yaw prediction error curve下载:
全尺寸图片
图5为航向模型参数辨识结果预报航向角速度曲线和实验测量角速度曲线。图6中给出了参数辨识结果角速度预报误差。从0~600 s预报误差曲线对比可知，FRDLS、MILS、FFLS相比于LS、LSSVM算法预报误差下降更快，有更快的辨识收敛速度。为了定量描述参数辨识精度，采用状态预报均方根误差(root mean square error，RMSE)作为航向模型辨识结果的精度指标，可根据式(23)求航向角速度预报均方根误差。
$$ {\text{RMSE}} = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {\boldsymbol e_i^2} } $$ (23) 式中
${\boldsymbol e_i}$ 为式(21)中第i组数据航向角速度预报误差，角速度预报误差结果和辨识过程收敛性见表2。表 2 角速度预报误差Table 2 Yaw prediction error算法 RMSE/(°·s−1) 收敛性 LS 0.043 9 收敛 FFLS 0.037 4 发散 MILS 0.022 1 发散 LSSVM — 发散 FRDLS 0.021 5 收敛 根据表2，FFLS相比于LS算法，引入遗忘因子后参数辨识精度有所提高，同时导致辨识过程产生了发散现象；FRDLS算法引入满秩分解后，提升了FFLS算法辨识精度的同时，消除了辨识过程中的发散，说明遗忘因子和满秩分解两种新息处理方式相结合具有合理性。LSSVM算法在1800 s后发散部分较多，对应的预报均方根误差结果为94943°·s−1，FRDLS算法在1800 s后未出现参数发散现象，说明了FRDLS算法在数据欠激励、环境扰动工况下具有更好的收敛性。
3.3 算法实时性分析
辨识船舶航向模型参数的过程中，使用了179 996组测量数据对参数估计结果进行了更新，采用Qt软件C++编程语言，在window10系统下记录了FFLS、MILS、LSSVM和 FRDLS算法每次使用测量数据从计算数据向量到参数辨识值更新的时间，得到算法单次循环耗时曲线如图7所示。
图 7 算法单次循环耗时曲线Fig. 7 Algorithm time consumption for a single loop下载:
全尺寸图片
算法单次循环最大耗时和平均耗时见表3。从算法耗时统计结果可知，FRDLS算法单次循环平均耗时为27.39
$ {\text{μ s}} $ ，远小于LSSVM算法，耗时介于FFLS算法和MILS算法之间。由于满秩分解过程会额外消耗算法单次循环时间，牺牲了FRDLS算法部分实时性，以提升收敛性，耗时为FFLS算法的3倍左右，实时性在可接受范围以内。表 3 算法单次循环耗时Table 3 Algorithm time consumption for a single loop算法 最大耗时 $/{\text{ms} }$ 平均耗时 $/{\text{μs} }$ FFLS 1.149 8.304 MILS 10.33 57.44 LSSVM 103.1 572.5 FRDLS 4.480 27.39 4. 结束语
本文针对数据欠激励受扰动系统参数在线辨识算法实时性、参数收敛性和辨识精度难以同时保证的问题，提出了一种满秩分解最小二乘算法，并使用该算法对长时间数据欠激励和海洋环境扰动工况下的船舶航向运动模型参数进行辨识。辨识结果表明，存在数据欠激励时，该算法能抑参数发散，并获得更高的辨识精度。且改进后的算法仍然计算量小、实现简单、实时性强、通用性好，具有较高实用价值，为实现准确稳定的航向自适应控制和欠驱动船自动靠泊奠定了基础。
-
图 1 实验数据曲线
Fig. 1 Experiment data curves
下载:
全尺寸图片
图 2 K辨识结果曲线
Fig. 2 K-coefficient curves
下载:
全尺寸图片
图 3 T辨识结果曲线
Fig. 3 T-coefficient curves
下载:
全尺寸图片
图 4
${\delta _d} $ 辨识结果曲线Fig. 4
${\delta _d}$ -coefficient curves下载:
全尺寸图片
图 5 角速度预报曲线
Fig. 5 Yaw prediction curve
下载:
全尺寸图片
图 6 角速度预报误差曲线
Fig. 6 Yaw prediction error curve
下载:
全尺寸图片
图 7 算法单次循环耗时曲线
Fig. 7 Algorithm time consumption for a single loop
下载:
全尺寸图片
表 1 参数辨识结果
Table 1 Parameters identification results
算法 $ K $ $T/{\text{s} }$ ${\delta _d}$/(°) LS 0.012 59.5 −0.23 FFLS 0.051 93.2 −0.31 MILS 0.059 116.8 −0.20 LSSVM 0.055 94.3 −0.33 FRDLS 0.038 70.9 −0.36 表 2 角速度预报误差
Table 2 Yaw prediction error
算法 RMSE/(°·s−1) 收敛性 LS 0.043 9 收敛 FFLS 0.037 4 发散 MILS 0.022 1 发散 LSSVM — 发散 FRDLS 0.021 5 收敛 表 3 算法单次循环耗时
Table 3 Algorithm time consumption for a single loop
算法 最大耗时 $/{\text{ms} }$ 平均耗时 $/{\text{μs} }$ FFLS 1.149 8.304 MILS 10.33 57.44 LSSVM 103.1 572.5 FRDLS 4.480 27.39 -
[1] 张笛, 赵银祥, 崔一帆, 等. 智能船舶的研究现状可视化分析与发展趋势[J]. 交通信息与安全, 2021, 1(29): 7–16 ZHANG Di, Zhao Yinxiang, CUI Yifan, et al. A Visualization Analysis and Development Trend of Intelligent Ship Studies[J]. Journal of transport information and safety, 2021, 1(29): 7–16. [2] 封波. 智能船舶发展战略规划研究[J]. 船舶工程, 2020, 42(3): 1–8. https://www.cnki.com.cn/Article/CJFDTOTAL-CANB202003011.htm FENG Bo. Research on intelligent ship development strategic planning[J]. Ship engineering, 2020, 42(3): 1–8. https://www.cnki.com.cn/Article/CJFDTOTAL-CANB202003011.htm [3] 严新平, 刘佳仑, 范爱龙, 等. 智能船舶技术发展与趋势简述[J]. 船舶工程, 2020, 42(3): 15–20. https://www.cnki.com.cn/Article/CJFDTOTAL-CANB202003008.htm YAN Xinping, LIU Jialun, FAN Ailong, et al. The development and tendency of intelligent vessel techniques[J]. Ship engineering, 2020, 42(3): 15–20. https://www.cnki.com.cn/Article/CJFDTOTAL-CANB202003008.htm [4] 徐海祥, 朱梦飞, 余文曌, 等. 面向智能船舶的自动靠泊鲁棒自适应控制[J]. 华中科技大学学报(自然科学版), 2020, 48(3): 25–40. https://www.cnki.com.cn/Article/CJFDTOTAL-HZLG202003005.htm XU Haixiang, ZHU Mengfei, YU Wenzhao, et al. Robust adaptive control for automatic berthing of intelligent ships[J]. Journal of Huazhong University of Science and Technology(natural science edition), 2020, 48(3): 25–40. https://www.cnki.com.cn/Article/CJFDTOTAL-HZLG202003005.htm [5] PAN Wei, XIE Xinlian, HE Ping, et al. An automatic route design algorithm for intelligent ships based on a novel environment modeling method[J]. Ocean engineering, 2021, 237:109603. [6] 秦贝贝, 陈增强, 孙明玮, 等. 基于自适应神经模糊推理系统的船舶航向自抗扰控制[J]. 智能系统学报, 2020, 15(2): 255–263. https://www.cnki.com.cn/Article/CJFDTOTAL-ZNXT202002009.htm QIN Beiei, CHEN Zengqiang, SUN Mingwei, et al. Active disturbance rejection control of ship course based on adaptive-network-based fuzzy inference system[J]. CAAI transactions on intelligent systems, 2020, 15(2): 255–263. https://www.cnki.com.cn/Article/CJFDTOTAL-ZNXT202002009.htm [7] 杨迪, 郭晨, 朱玉华, 等. 欠驱动船舶神经网络自适应路径跟踪控制[J]. 智能系统学报, 2018, 13(2): 254−260 YANG Di, GUO Chen, ZHU Yuhua, et al. Neural network adaptive path tracking control for underactuated ships[J]. CAAI transactions on intelligent systems, 2018, 13(2): 254−260. [8] LIU Zhixiang, ZHANG Youmin, YU Xiang, et al. Unmanned surface vehicles: an overview of developments and challenges[J]. Annual reviews in control, 2016, 41: 71–93. doi: 10.1016/j.arcontrol.2016.04.018 [9] FENG Ding, LIU P X, LIU Guangjun. Multi-innovation least-squares identification for system modeling[J]. IEEE transactions on systems, man, and cybernetics, part B (cybernetics), 2010, 40(3): 767–778. doi: 10.1109/TSMCB.2009.2028871 [10] 孙明轩, 毕宏博. 学习辨识: 最小二乘算法及其重复一致性[J]. 自动化学报, 2012, 38(5): 698–706. doi: 10.3724/SP.J.1004.2012.00698 SUN Mingxuan, BI Hongbo. Learning identification: least squares algorithms and their repetitive consistency[J]. Acta automatica sinica, 2012, 38(5): 698–706. doi: 10.3724/SP.J.1004.2012.00698 [11] 秦余钢, 马勇, 张亮, 等. 基于改进最小二乘算法的船舶操纵性参数辨识[J]. 吉林大学学报(工学版), 2016, 46(3): 897–903. https://www.cnki.com.cn/Article/CJFDTOTAL-JLGY201603033.htm QIN Yugang, MA Yong, ZHANG Liang, et al. Parameter identification of ship’s maneuvering motion based on improved least square method[J]. Journal of Jilin University(engineering and technology edition), 2016, 46(3): 897–903. https://www.cnki.com.cn/Article/CJFDTOTAL-JLGY201603033.htm [12] SONG Chunyu, ZHANG Xianku, ZHANG Guoqing. Nonlinear identification for 4 DOF ship maneuvering modeling via full-scale trial data [J]. IEEE transactions on industrial electronics, 2021, 99:1. [13] 孙功武, 谢基榕, 王俊轩. 基于动态遗忘因子递推最小二乘算法的船舶航向模型辨识[J]. 计算机应用, 2018, 38(3): 900–904. doi: 10.11772/j.issn.1001-9081.2017082041 SUN Gongwu, XIE Jirong, WANG Junxuan. Ship course identification model based on recursive least squares algorithm with dynamic forgetting factor[J]. Journal of computer applications, 2018, 38(3): 900–904. doi: 10.11772/j.issn.1001-9081.2017082041 [14] ZHAO Baigang, ZHANG Xianku. An improved nonlinear innovation-based parameter identification algorithm for ship models[J]. The journal of navigation, 2021, 74(3): 549–557. doi: 10.1017/S0373463321000102 [15] 史颂辉, 丁世飞. 基于能量的结构化最小二乘孪生支持向量机[J]. 智能系统学报, 2020, 13(5): 1013–1019. https://www.cnki.com.cn/Article/CJFDTOTAL-ZNXT202005027.htm SHI Songhui, DING Shifei. Energy-based structural least square twin support vector machine[J]. CAAI transactions on intelligent systems, 2020, 13(5): 1013–1019. https://www.cnki.com.cn/Article/CJFDTOTAL-ZNXT202005027.htm [16] CHEN Sugen, WU Xiaojun. Multiple birth least squares support vector machine for multi-class classification[J]. International journal of machine learning and cybernetics, 2017, 8(6): 1731–1742. doi: 10.1007/s13042-016-0554-7 [17] LAURAIN V, TOTH R, PIGA D, et al. An instrumental least squares support vector machine for nonlinear system identification[J]. Automatica, 2015, 54: 340–347. doi: 10.1016/j.automatica.2015.02.017 [18] WANG Zihao, ZOU Zaojian, SOARES C G. Identification of ship manoeuvring motion based on nu-support vector machine[J]. Ocean engineering, 2019, 183: 270–281. doi: 10.1016/j.oceaneng.2019.04.085 [19] 谢朔, 初秀明, 柳晨光, 等. 基于改进 LSSVM 的船舶操纵运动模型在线参数辨识方法[J]. 中国造船, 2018, 59(2): 178–188. doi: 10.3969/j.issn.1000-4882.2018.02.019 XIE Su, CHU Xiuming, LIU Chenguang, et al. Online parameter identification method for ship maneuvering models based on improved LSSVM[J]. Ship building of China, 2018, 59(2): 178–188. doi: 10.3969/j.issn.1000-4882.2018.02.019 [20] ZHU Man, SUN Wuqiang, HAHN AXEL, et al. Adaptive modeling of maritime autonomous surface ships with uncertainty using a weighted LS-SVR robust to outliers[J]. Ocean engineering, 2020, 200:107053. [21] LUO Weilin, MOREIRA L, SOARES C G. Manoeuvring simulation of catamaran by using implicit models based on support vector machines[J]. Ocean engineering, 2014, 82: 150–159. doi: 10.1016/j.oceaneng.2014.03.008 [22] LUO Weilin, MOREIRA L , SOARES C G . Modeling of ship maneuvering motion using neural networks[J]. Journal of marine science and application, 2016, 15(4): 426–432. doi: 10.1007/s11804-016-1380-8 [23] RAJESH G, BHATTACHARYYA S K. System identification for nonlinear maneuvering of large tankers using artificial neural network[J]. Applied ocean research, 2008, 30(4): 256–263. doi: 10.1016/j.apor.2008.10.003 [24] WANG Ning, MENG JOO Er , HAN Min. Large tanker motion model identification using generalized ellipsoidal basis function-based fuzzy neural networks[J]. IEEE transactions on cybernetics, 2015, 45(12): 2732–2743. doi: 10.1109/TCYB.2014.2382679 [25] 褚式新, 茅云生, 董早鹏, 等. 基于极大似然法的高速无人艇操纵响应模型参数辨识[J]. 兵工学报, 2020, 41(1): 127–134. https://www.cnki.com.cn/Article/CJFDTOTAL-BIGO202001015.htm CHU Shixin, MAO Yunsheng, DONG Zaopeng, et al. Parameter identification of high-speed USV maneuvering response model based on maximum likelihood algorithm[J]. Acta armamentar II, 2020, 41(1): 127–134. https://www.cnki.com.cn/Article/CJFDTOTAL-BIGO202001015.htm [26] 张炜灵, 蔡烽, 王骁. 基于解析模型的小样本操纵性KT指数辨识方法[J]. 中国航海, 2020, 43(3): 62–67. doi: 10.3969/j.issn.1000-4653.2020.03.011 ZHANG Weiling, CAI Feng, WANG Xiao. Analytical model for maneuverability KT index identification with small sample[J]. Navigation of China, 2020, 43(3): 62–67. doi: 10.3969/j.issn.1000-4653.2020.03.011 [27] 徐健平. 船舶自适应舵的一种算法研究[D]. 大连: 大连海事大学, 2012.XU Jianping. The research of a algorithm of ship adaptive autopilot [D]. Dalian: Dalian Maritime University, 2012.
