CN104536944B - 基于改进的极大似然方法和置信压缩滤波的船舶参数辨识方法 - Google Patents

Abstract

本发明公开了基于改进的极大似然方法和置信压缩滤波的船舶参数辨识方法。包括以下几个步骤：步骤一，测量执行机构对船舶施加的力和力矩作为控制输入，采集船舶的位移和姿态角信息作为输出观测量，设定初始状态估计值，根据控制输入和输出观测量建立船舶模型，得到船舶模型中待辨识的参数；步骤二，根据船舶当前的状态估计值和输出观测量η，利用改进的极大似然方法，辨识船舶模型的参数估计值；步骤三，根据辨识出的船舶模型的参数估计值，利用置信压缩滤波器，更新当前的状态估计值；重复步骤二～步骤三，直到停止输出船舶的输出观测量，得到每个采样时刻的参数估计值。本发明计算简单，执行效率高，辨识结果准确，能够在线运行。

Description

基于改进的极大似然方法和置信压缩滤波的船舶参数辨识 方法

技术领域

本发明属于一种船舶参数辨识方法，尤其涉及一种针对船舶时域状态空间模型中未知参数的估计的，基于改进的极大似然方法和置信压缩滤波的船舶参数辨识方法。

背景技术

在航行过程中船舶的姿态控制是一个复杂的问题。船舶有6个自由度，船舶的运动方程存在高度的非线性，水动力系数之间相互耦合；同时，船舶在航行过程中受到不确定的洋流、海浪、风等的干扰，执行机构存在滞后性，这样就对船舶的控制造成了很大困难。为了对船舶进行准确的控制，一般需要建立准确的船舶模型，进一步需要对包括水动力系数在内的船舶参数进行辨识，减小控制的不确定性。船舶参数可以通过实测或理论估计得到，前者往往要借助于特殊的仪器，成本较高，而且时间消耗较大，获取麻烦；后者采用经验或切片法等理论，往往不够准确。通过辨识方法确定船舶参数综合了二者的优点，有较高的研究价值。参数辨识还可以用于解决其它领域的一些测试方法中要求的物理量难以直接获取的情况。

极大似然方法是一种迭代算法，能够渐进地取得极大值及其对应的参数值，这样就能够辨识得到船舶参数。然而极大似然方法需要经过相当多步迭代后才能收敛到正确结果，收敛速度慢。

发明内容

本发明的目的是提供一种能够在线运行的基于改进的极大似然方法和置信压缩滤波的船舶参数辨识方法。

本发明是通过以下技术方案实现的：

基于改进的极大似然方法和置信压缩滤波的船舶参数辨识方法，包括以下几个步骤：

步骤一：测量执行机构对船舶施加的力和力矩作为控制输入τ，采集船舶的位移和姿态角信息作为输出观测量η，设定初始状态估计值，根据控制输入τ和输出观测量η建立船舶模型，得到船舶模型中待辨识的参数θ；

步骤二：根据船舶当前的状态估计值和输出观测量η，利用改进的极大似然方法，辨识船舶模型的参数估计值；

步骤三：根据辨识出的船舶模型的参数估计值，利用置信压缩滤波器，更新当前的状态估计值；

步骤四：重复步骤二～步骤三，直到停止输出船舶的输出观测量，得到每个采样时刻的参数估计值。

本发明基于改进的极大似然方法和置信压缩滤波的船舶参数辨识方法还包括：

1、船舶模型中待辨识的参数θ为：

其中，待辨识的参数θ包括的各项通过A_x表示A对下标x方向上的附加质量，A＝X, Y,N定义

2、根据船舶当前的状态估计值及输出观测量η，利用改进的极大似然方法，辨识船舶模型的参数估计值的方法为：

(1)令k＝1，初始化θ₀；

(2)根据当前时刻k的状态估计值输出观测值η_k，利用船舶模型求得概率其中θ_k-1表示k-1时刻得到的参数值；

(3)计算概率i＝1,2,...,N，N为所用参数值θ_k的粒子数，θ_k表示k时刻得到 的参数值，定义概率

利用先验知识更新θ_k，建立一阶自回归模型：

θ_k＝θ_k-1+ε；

其中ε服从方差收敛的高斯分布或者ε服从边界收敛的均匀分布，

在时刻k，根据一阶自回归模型生成一组支持点对于每个计算得到概率

(4)计算概率的最大值得到极大似然概率

其中，为k时刻得到的参数θ_k的最优的估计值

(5)令k＝k+1，并从第(2)步迭代计算，直至k＝M+1，M为输出观测值的个数。

3、根据辨识出的船舶模型的参数估计值，利用置信压缩滤波器，更新当前的状态估计值的方 法为：

(1)首先选取离散支持点数m_k；

(2)求取支持点s_k,i,i＝1,2,...m_k，支持点s_k,i满足等式：

其中p(η_k|ξ_k),p(ξ_k/ξ_k-1＝s_k-1,i)；

(3)求取当前的状态估计值为：

有益效果：

本发明采用了改进的极大似然方法，无需对观测数据进行批处理，无需进行大量迭代，在每个采样时刻都能取得一个以极大似然概率保证的参数估计值；采用了置信压缩滤波算法进行参数辨识，能够用于复杂的非线性状态空间模型，具有很高的滤波精度，避免了粒子滤波这样的滤波方法引入的粒子退化、多样性匮乏等问题；计算简单，执行效率高，辨识结果准确，能够在线运行。

附图说明

图1是整个船舶参数辨识方法原理图；

图2是艏摇角速度r的仿真结果；

图3是和艏摇角ψ的仿真结果；

图4是参数辨识结果；

图5是置信压缩滤波结果。

具体实施方式

下面将结合附图对本发明做进一步详细说明。

本发明描述的是一种新的船舶参数辨识方法，利用了改进的极大似然方法和置信压缩滤波方法，方案图如图1所示。具体步骤如下：

(1)利用船舶自身的推进器、舵等执行机构对船舶施加力和力矩，使船舶做特定运动，具体为：

通过控制器对船舶的螺旋桨、舵等执行机构进行控制，记录下每个采样时刻的推进力及力矩；船舶的整体运动方式应该严格按照船舶模型确定；为了使船舶完成特定的运动方式，可以采用额外的机械结构对船舶施加外部力和力矩，并且所施加的外部力及力矩均要求可测 并且进行记录。

(2)利用测量系统测出地坐标系下，船舶在每个离散采样时刻的位移和姿态角信息，其为船舶模型的输出观测值，具体为：

这里所说的船舶模型为：

根据牛顿定理和拉格朗日公式，船体的6自由度非线性动态方程可以描述为

以上动态方程同时适用于其它自由度。其中，M＝M_RB+M_A为惯性矩阵，M_RB是船体的刚体质量矩阵，M_A是附加质量矩阵；C(ξ)＝C_RB(ξ)+C_A(ξ)为由地球自转引起的科氏力和向心力矩阵，C_RB(ξ)是刚体科氏力和向心力矩阵，C_A(ξ)是附加质量科氏力和向心力矩阵；D(ξ)为流体阻尼矩阵，包括船体震动引起的阻尼、摩擦阻尼、洋流阻尼、涡流阻尼等；g(η)为重力及浮力同时作用引起的恢复力及力矩；τ为控制输入，这里忽略了环境干扰力等外部作用。观测变量η表示船体在地坐标系中的位置与姿态，用体坐标系坐标原点在地坐标系中的坐标(x,y,z)表示船体在地坐标系中的位置，同时用地坐标系与体坐标系之间的欧拉角(φ,θ,ψ)表示船体的姿态；状态变量ξ表示船体的体坐标系线速度和角速度，其中(u,v,w)表示船体的线速度，(p,q,r)表示船体绕各轴旋转的角速度。

通常船体的6自由度模型非常复杂，在其基础上忽略垂荡、纵摇、横滚运动，即认为w＝p＝q＝0，同时考虑水下船体的工作情况和对称性，建立水下船体的平面3自由度模型。假设船体具有均匀质量分布，关于xz平面对称，则I_xy＝I_yz＝0；假设浮心和重心重合，这样浮力和重力作用于一点。

根据式(1)和式(2)，重新定义ξ＝[u v r^T]，η＝[x y ψ]^T。M_RB是一个对称阵， m为船体质量，包含自由流动空间的水质量，x_G和y_G为水平面 上的船体的重心位置。矩阵中表示A 对下标x方向上的附加质量，定义为例如其它类似。

矩阵C(ξ)是一个反对称阵，其中 矩阵矩阵中带 下标的元素与矩阵M_A中元素的意义相同，例如

向量g(η)表示重力和浮力及相应力矩，这里认为船体能够稳定在某一水深位置，取g(η)＝[0 0 0]^T。

向量τ表示外力及外力矩作用，不包括环境干扰，τ＝[τ_X τ_Y τ_Z]^T。为了评估参数辨识的效果，这里认为所有施加的外力及外力矩均是可测的。

欧拉变换方程将体坐标系中观察到的物理量：前进线速度、侧移线速度、艏摇角速度变换到地坐标系中，得到x轴上速度、y轴上速度及艏摇角速度

参数：

将式(1)和式(2)进一步表示为

设采样时间间隔为Δ，采用一阶欧拉法进行离散化，得到船舶的非线性离散模型

ξ_k＝ξ_k-1+ΔM^-1[τ_k-1-C(ξ_k-1)ξ_k-1-D(ξ_k-1)ξ_k-1-g(η_k-1)] (5)

η_k＝η_k-1+ΔJ(η_k-1)ξ_k (6)

船舶模型的输出观测值指每个时刻的η_k。

(3)利用改进的极大似然方法，根据船舶的状态估计值及输出观测值辨识船舶得到每个 时刻的船舶参数估计值，具体为：

这里的所说的改进的极大似然方法，指的是：

将式(5)和式(6)写成含参非线性系统通用形式如下

其中，{ξ_k}是一阶马尔科夫过程；η_k的条件概率分布仅仅与ξ_k有关；θ为未知参数；w_k-1和v_k分别为过程噪声和输出噪声，且统计特性已知，这里假设w_k-1和v_k均服从均值为0的高斯分布，即w_k-1～N(0,W)，v_k～N(0,V)。已知输出方程得到的η_k的观测序列η_1:M＝{η_k,k＝1,2,...,M}，其中M为总的观测值的数量。这里需要解决两个问题，即参数θ的辨识问题和k时刻状态ξ_k的估计问题。

构造极大似然概率函数，非线性参数辨识问题描述为参数θ下η_1:M的极大似然条件概率问题，即

其中p_θ(η_1:M)＝p(η_1:M|θ)，直接计算其最大值对于非线性函数来说，往往是不能实现的。利用贝叶斯公式，可以得到

式(9)认为η_k之间相互独立，其表明，在每个时刻k都取得极大值是取得极大值 的充要条件。为了求得只需得到k时刻的极大值当所有M个观测时刻都取得似然函数极大值时，就可以取得p_θ(η_1:M)的极大值，进而取得参数θ的估计值。在执行EM算法时，只需在当前时刻进行一步迭代，利用M个观测值进行共M次迭代即可。所提出的算法因为每个时刻只有一次迭代，无须设置迭代结束条件。具体来说，在一次极大 似然计算时，只需利用k时刻的观测值η_k，计算出当前的似然概率即可，所研究的问题等价转化为

利用Chapman-Kolmogorov方程，对p(η_k|θ)给出数值解法

p(η_k|θ)＝∫p(η_k|ξ_k)p(ξ_k|θ)dξ_k (11)

直接对式(11)进行积分计算仍然不易实现。进一步利用蒙特卡罗方法，从p(ξ_k|θ)中抽取 可以得到

式(12)由大数定理保证成立。由于p(ξ₀|θ)的先验分布是已知的，根据式(12)就可以得到p(η₁|θ)。为了完成参数θ的迭代计算，必须对下一时刻的p(ξ₁|θ)进行估计。由于式(7)为非线性状态方程，这里利用置信压缩滤波计算式(11)中p(ξ_k|θ)(k＝1,2,...,M)的先验概率分布。

定义在含M个观测值的情况下，只要满足似然概率就将得到的θ_k作为第k步迭代的参数估计值，其中θ_k-1表示k-1时刻得到的参数值，θ_k表示k时刻得到的参数新值。具体算法如下：

1)令k＝1，初始化θ₀

2)期望阶段：

计算

计算N为所用θ_k的粒子数

其中按式(13)～(15)产生。

3)最大化阶段：

极大似然概率为得到相应 的参数估计值

4)令k＝k+1，并从第(2)步迭代计算，直至k＝M+1

下面给出期望和最大化两阶段算法的具体实现。

1计算期望

注意到式(12)并不是关于θ的显式表达式，直接求梯度进行最大期望计算是不能实现的。设参数存在于一定范围内，θ∈[-B,B],B∈R⁺。对于时刻k，生成一组支持点在这N个支持点确定的可行解中找到最优值，将极值问题转化为有界集合内的参数寻优问题。

利用先验知识可以更准确地估计θ_k，为此考虑利用前一步参数的极大似然估计结果θ_k-1。下面给出在θ_k-1的基础上更新θ_k的计算公式。考虑如下一阶自回归模型

θ_k＝θ_k-1+ε (13)

1)ε服从方差收敛的高斯分布

ε～N(0,b_k) (14)

2)ε服从边界收敛的均匀分布

ε～U[-b_k,b_k] (15)

两种情况下，θ_k均收敛，其中b_k＝k^-αb_k-1，b₀＝b，α(0＜α≤1)为衰减因子，b∈R⁺。

在时刻k，根据式(13)和以上模型生成一组支持点对于每个按照式(12)计算得到则取其中的最大值，即

这样即可完成θ_k的收敛的迭代更新和极大似然函数的求取。

2求取最大化的参数值

满足式(16)的为k时刻得到的参数θ_k的最优的估计值即

改进的极大似然方法的思想是利用期望-最大化算法，使用搜索算法，在参数的可行解空间中确定出参数值。每个参数产生N个支持点，用于单参数辨识时，搜索空间是1×N维的；用于n个参数辨识时，搜索空间是Nⁿ维。

利用改进的极大似然方法，根据当前时刻k的状态估计值输出观测值η_k辨识得到当前 时刻船舶参数的估计值

(4)利用置信压缩滤波器，根据每个时刻的受力信息、船舶的输出观测值和前一时刻辨识的参数值，估计船舶的状态值，具体为：

这里所说的置信压缩滤波器，指的是：

对于式(5)式(6)表示的非线性系统，由步骤(3)辨识出参数前一时刻值后，可以表示成下式

其中，

滤波就是根据η_1:k得到ξ_k的估计值。如果得到概率密度函数p(ξ_k|η_1:k)，ξ_k的估计值将很容易得到。设p(ξ₀|η₀)＝p(ξ₀)已知，即ξ₀的先验分布已知。

定义置信压缩滤波的离散实现中似然概率分布F(ξ)∈P，是对真实分布F的估计；定义映射φ，其将分布由p_k-1＝p(ξ_k-1|η_1:k-1)映射到p_k＝p(ξ_k|η_1:k)，即p_k＝φ(p_k-1)。置信压缩滤波离散实现的算法为：

1)初始化

选取使之服从ξ₀的先验分布

2)状态估计

For k＝1,2,...M

1.选取m_k作为离散支持点数

2.计算离散支持点s_k,i(i＝1,2,…m_k)：

For i＝1,2,...m_k，按照下式计算得到s_k,i

其中

3.按照下式近似

是指示函数，定义为

利用船舶的受力τ、输出值η及待辨识参数的前一时刻的估计值，由置信压缩滤波估计出当前时刻k的状态估计值

(5)在观测数据结束前，不断执行(3)和(4)两步，具体为：

只要还有船舶输出观测数据，就不断执行(3)和(4)两步的计算，即结束条件为k＝M。算法的整个过程的计算是收敛的、稳定的，参数估计值将不断趋于正确值。

实验

通过一个水下无人器ODIN进行实验，ODIN(Omni-directional IntelligentNavigator,全方位智能导航器)是由Hawaii大学研制的具有6个自由度的无人水下航行器。

该潜器具有对称性，即船体具有均匀质量分布。

通过推进器，使ODIN在水平面上作绕垂线的旋转运动，即可忽略垂荡、纵摇、横滚运动。

同时考虑由于流体动力学引入的附加质量，考虑海浪干扰影响。在上述条件下，推导出ODIN的离散化模型。

ψ_k＝ψ_k-1+Δr_k+n_o,k (23)

式(22)和式(23)确定的离散模型中，N_x(x＝r,r|r|)表示N对下标x方向上的附加 质量，定义为

其中：N____关于z轴旋转(艏摇)力矩，r___绕z轴旋转的(艏摇)角速度。

N_r及N_r|r|为待辨识参数。N_r和N_r|r|的真实值分别为-30和-80。

取r和ψ的初始值为零，Δ＝0.05s，n_s,k-1～N(0,0.0001)，n_o,k～N(0,0.01)，固定力矩 τ_p＝30N，进行5s仿真，得到101个采样时刻。艏摇角速度r和艏摇角ψ的仿真结果分别如图2和图3所示。

利用改进的极大似然方法，对水下机器人进行参数辨识及滤波。取初始状态方差P₀＝0.1，N_r和N_r|r|的初值分别设为-45和-60，期望-最大化算法中取α＝0.015，利用方差收敛的高斯分布，其初始方差设为5。参数辨识和滤波所用支持点均取100。参数辨识结果如图4所示，置信压缩滤波结果如图5所示。从图4可以看出，参数辨识过程收敛，N_r和N_r|r|的参数辨识结果分别为-26.6和-79.6，接近真实值。

Claims (2)

1.基于改进的极大似然方法和置信压缩滤波的船舶参数辨识方法，其特征在于，包括以下几个步骤：

步骤一：测量执行机构对船舶施加的力和力矩作为控制输入τ，采集船舶的位移和姿态角信息作为输出观测量η，设定初始状态估计值，根据控制输入τ和输出观测量η建立船舶模型，得到船舶模型中待辨识的参数θ；

步骤二：根据船舶当前的状态估计值和输出观测量η，利用改进的极大似然方法，辨识船舶模型的参数估计值；

步骤三：根据辨识出的船舶模型的参数估计值，利用置信压缩滤波器，更新当前的状态估计值；

步骤四：重复步骤二～步骤三，直到停止输出船舶的输出观测量，得到每个采样时刻的参数估计值；

所述的船舶模型中待辨识的参数θ为

其中，待辨识的参数θ包括的各项通过A_x表示A对下标x方向上的附加质量，A＝X,Y,N定义

所述的根据船舶当前的状态估计值及输出观测量η，利用改进的极大似然方法，辨识船舶模型的参数估计值的方法为：

(1)令k＝1，初始化θ₀；

(2)根据当前时刻k的状态估计值输出观测值η_k，利用船舶模型求得概率其中θ_k-1表示k-1时刻得到的参数值；

(3)计算概率N为所用参数值θ_k的粒子数，θ_k表示k时刻得到的参数值，定义概率

利用先验知识更新θ_k，建立一阶自回归模型：

θ_k＝θ_k-1+ε；

其中ε服从方差收敛的高斯分布或者ε服从边界收敛的均匀分布，

在时刻k，根据一阶自回归模型生成一组支持点对于每个计算得到概率

(4)计算概率的最大值得到极大似然概率

其中，为k时刻得到的参数θ_k的最优的估计值

(5)令k＝k+1，并从第(2)步迭代计算，直至k＝M+1，M为输出观测值的个数。

2.根据权利要求1所述的极大似然方法和置信压缩滤波的船舶参数辨识方法，其特征在于：所述的根据辨识出的船舶模型的参数估计值，利用置信压缩滤波器，更新当前的状态估计值的方法为：

(1)首先选取离散支持点数m_k；

(2)求取支持点s_k,i,i＝1,2,...m_k，支持点s_k,i满足等式：

其中p(η_k|ξ_k),p(ξ_k/ξ_k-1＝s_k-1,i)；

(3)求取当前的状态估计值为：

是指示函数，定义为