一种机动浮标在随机海浪干扰下的自守卫控制方法与流程

本发明属于海洋浮标领域，具体涉及一种机动浮标在随机海浪干扰下的自守卫控制方法。

背景技术：

海洋浮标是以锚定在海上的观测浮标为主体组成的海洋水文水质气象自动观测站。它能按规定要求长期、连续地收集所需海洋水文水质气象资料，特别是能收集到调查船难以恶劣天气及海况的资料。机动浮标是一种特殊的海洋浮标，在受到外界环境影响而漂离到给定的守卫区域外时，它可以依靠自身的动力，自主航行到给定的守卫点，继续执行自守卫测量任务。但是，当机动浮标在近水面自守卫航行时，在随机海浪的干扰下，由于海浪与机动浮标运动之间的相互作用产生的波浪增阻会降低机动浮标的主机推进功率，使得机动浮标在自守卫航行过程中的能量严重损耗，进而降低机动浮标的航行效率。

技术实现要素：

本发明的目的在于提供一种机动浮标在随机海浪干扰下的自守卫控制方法，其能有效减少机动浮标在自守卫航行时的波浪增阻，进而减少能量损耗，提高机动浮标的航行效率。

一种机动浮标在随机海浪干扰下的自守卫控制方法，具体包括以下步骤：

步骤1、根据机动浮标在随机海浪干扰下的纵摇、横摇、垂荡运动的时间序列，采用mar(multiautoregressivemodel)模型法计算机动浮标在上述三自由度上运动的交叉谱；

步骤2、根据机动浮标自身的模型参数通过matlab中的mss工具包计算出浮标在上述三自由度上运动的rao(responseamplitudeoperator)数据；

步骤3、根据步骤1中获得的交叉谱和步骤2中获得的rao数据计算得到机动浮标的海浪方向谱，三者之间的关系为：

式中，φij(ω)为通过机动浮标的时间序列获取的交叉谱，raoi(ω,θ)表示在第i个自由度上，波浪入射角为θ，频率为ω时的响应幅值算子，s(ω,θ)为海浪方向谱；

步骤4、根据机动浮标的实际航行情况，在给定的自守卫区域内，设置虚拟节点，每个虚拟节点中都包含此时刻的浪向角信息；

步骤5、根据机动浮标在随机海浪干扰下受到的波浪增阻与浪向角之间的函数关系，将机动浮标在自守卫航行时的能量损耗作为目标函数；

步骤6、采用a*算法对机动浮标在随机海浪干扰下自守卫航行时，从守卫区域边缘点到达给定守卫点进行路径规划，使得机动浮标在规划的航路上航行时能量损耗最低。

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，步骤1中mar模型建立方法为：

步骤1、假设{zk}是k维随机平稳的零均值时间序列，表示如下：

zk＝[z1z2...zk]^t

其中，zi(s·δt),s＝1,2,...,n，i＝1,2,...,k,其中δt为采样时间，令采样时间δt＝1，则zi(s·δt)＝zi(s)；

步骤2、由于所需要分析的是机动浮标三自由度的运动时间序列，在每个时间序列中，其前p个时间序列值对第s个时间序列进行建模，则有：

z(s)-a1z(s-1)-...-ap(s-p)＝ε(s)

其中，ai为k×k维的权重系数矩阵；ε(s)＝[ε1ε2...εk]^t表示均值为0，方差为σi的高斯白噪声序列，由于εi与z(s)的取值无关，因此ε(s)与ε(s-p)的协方差为零；由此可得白噪声向量的协方差矩阵如下：

式中，σij表示εi与εj之间的协方差；将上式改写如下：

令：v(s)＝[z(s-1)z(s-2)...z(s-p)]^t

λp＝[λ1λ2...λp](k×(k·p))

则mar模型的最终形式为：

z(s)＝λpv(s)+ε(s)

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，对上述权重系数矩阵ai求解时采用whittle递推法；由于mar模型的权值系数aj(j＝1,2,...,p)随着模型阶数p的变化而改变；因此，令：λpj＝(ap1ap1...apj)；j＝1,2,...,p，ri(i＝1,2,...,p)为相关矩阵，对于不同迟滞l相关矩阵中的第i,j个元素的计算公式表示为：

因此，递推公式为：

针对上述递推法求出的p项权重系数矩阵采用最优阶数判断准则来选取最有效的一组权重系数矩阵；该准则的描述为：对于长度为n的k维平稳随机序列{zk}，其某一步预报误差方差矩阵的估计量用dp来表示，最终的预报误差为其行列式的值：

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，步骤1中计算机动浮标运动交叉谱的方法包括以下步骤：

根据上述建立的mar模型，选择适当的阶数，就可以准确地使用前p个历史值准确预测所求摇荡时序的当前值。

步骤1、由交叉谱密度的基本定义得到白噪声的交叉谱密度

式中f为变换后的频率，由于当且仅当s＝t时，即等式右边只是将σst进行傅里叶变换，而σst是一个常数项，对于常数项的傅里叶变换仍然是其本身；可得

步骤2、考虑残差的影响，令aij(0)＝-δij，其中：

根据上式可得：

其中，εi(s)是σj(s)的线性变换，模型当维数退化到一维时，aij(p)为脉冲响应函数：

此时，原始信号经过脉冲响应函数的线性变化后，在其g频点的功率谱密度phh(g)为：

phh(g)＝|b(g)|²pyy(g),-∞＜g＜∞

式中，pyy(g)是y(t)的功率谱密度，为脉冲响应函数h(τ)的频率响应函数；

分析可知，多维与一维的情况相同，但是多维误差的交叉谱为复数，可知：

其中，α,β＝1,2,...k；p(f)即为所求的交叉谱；考虑到其矩阵形式表示为：

π＝a(f)p(f)a(f)^h

式中，π为由σst组成的矩阵；假定矩阵a(f)为非奇异矩阵，上式两端同时乘以a(f)的逆矩阵和a(f)^h(即a(f)的共轭转置矩阵)，可得：

p(f)＝a(f)^-1πa(f)^-h

综上，便可获得机动浮标三自由度摇荡运动的交叉谱，最终可得海浪方向谱。

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，步骤4中虚拟节点的设置原则为：根据机动浮标在湖上和海上试验的实际情况，自守卫区域半径一般取50m左右，守卫内径在5m左右，而机动浮标自身的长度为1.35m。因此，本发明设置的虚拟节点是以自守卫内圆为中心点，在守卫区域内以5m的间距均匀设置虚拟节点。每个虚拟节点都包含此时的浪向角信息。

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，机动浮标在随机海浪干扰下产生的波浪增阻f(β)与浪向即遭遇角β之间的函数关系为：

f(β)＝p0+p1β+...+pkβi^k

式中，(p0,...,pk)表示所需要求出的多项式系数，在一级海况、二级海况、三级海况下，采用最小二乘法对试验数据进行拟合仿真，求得的系数矩阵分别为：p＝[2.73×10^-10-1.22×10^-71.85×10^-4-7.44×10^-4]^t，p＝[-1.24×10^-98.85×10^-7-2.22×10^-40.027]^t，p＝[-1.13×10^-91.53×10^-6-4.92×10^-40.07]^t。

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，步骤6具体包括以下步骤：

步骤6.1、设计代价函数，a*算法中的代价函数f(x)表示为：

f(x)＝g(x)+h(x)

式中，f(x)表示从起始节点n0开始到终止节点的最佳路径的代价总值；g(x)表示从起始节点n0开始到某一个节点x时的最佳路径代价值，h(x)表示从某一节点x到终止节点的实际最佳路径代价值，h(x)又被称为启发函数，当h(x)的值小于某一节点x到终止节点的实际距离时，可以得到最优解；

步骤6.2、根据波浪增阻与遭遇角之间的函数关系表达式，可得机动浮标航行至某一节点处的能量消耗值，设机动浮标的航速为v，起始节点的位置为(x0,y0)，在某一节点ni(xi,yi)节点处，航行至下一节点ni+1(xi+1,yi+1)的时间为t，则在这两个相邻节点的路径长度为l，则：

机动浮标在这段路径上航行的能量损失δploss(ni,ni+1)表示为：

把在整个规划的路径上损耗的总能量ploss作为代价函数，表示为：

步骤6.3、设计启发函数，选择欧几里得距离法作为机动浮标路径规划的启发方法，其对应的函数表达式表示为：

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，步骤6中a*算法具体包括以下步骤：

步骤1、建立三个表，分别为open表、closed表、parent表，其中，open表用来存放当前根据海流和海浪方向信息规划的节点信息、closed表用来存放已经选取的节点信息、parent表用来存放父节点信息；

步骤2、对算法进行初始化设置，即将open表、closed表和parent表清空，然后设置起始节点n0、终止节点nn的位置信息；

步骤3、将起始节点n0放入open表中，令g(n0)＝0，则f＝h；

步骤4、判断open表是否为空，若open表为空，则所需规划的路径不存在，算法结束，若open表不为空，则执行步骤5；

步骤5、把open表中的代价函数值最小的节点定义为nbest，判断nbest是否为终止节点nn，若nbest＝nn，则表示当前路径规划完成，算法结束，如果nbest≠nn，则执行步骤6；

步骤6、计算open表中所存入的节点的代价函数值ploss，并比较相对应的代价函数值的大小，把代价函数值最小的节点即nbest从open表中赋给closed表中；

步骤7、判断是否满足能量损耗最低的约束条件，若满足在该节点的能量损耗最低，则计算该节点的代价值，并与实际代价值进行比较，将代价值最小的节点添加到open表和parent表中，返回步骤6，若不满足能量损耗最低的约束条件，则返回步骤4。

本发明的有益效果在于：

本发明首先根据机动浮标在随机海浪干扰下的纵摇、横摇、垂荡运动的时间序列采用mar模型法获取机动浮标相对应运动的交叉谱φij(ω)，使用matlab获取机动浮标相应运动的rao数据，结合交叉谱和rao数据得到海浪方向谱s(ω,θ)，即可得到机动浮标自守卫航行时的实时浪向信息。然后根据机动浮标自守卫航行时的实际情况，合理设置虚拟节点，针对波浪增阻问题，通过设计a*算法对机动浮标自守卫航行进行能量最优航路规划。本发明可以有效减少机动浮标自守卫航行时的波浪增阻，进而减少能量损耗，提高航行效率。

附图说明

图1为本发明的机动浮标自守卫区域示意图；

图2为本发明的机动浮标的结构示意图；

图3为本发明的机动浮标在随机海浪干扰下的自守卫航行控制方法流程图；

图4为本发明的虚拟节点设置的示意图；

图5为本发明的a*算法的设计流程图。

具体实施方式

下面结合附图对本发明做进一步描述：

如附图1所示，为本发明的机动浮标自守卫区域示意图，如附图2所示，为本发明的机动浮标的结构示意图；机动浮标停机在守卫内圆中进行测量任务时，在外界未知环境干扰下，会逐渐漂流到守卫外圆之外。此时，机动浮标开启自守卫任务，即测量并计算从当前位置到达守卫点的目标航向与机动浮标艏向的差值，启动推进装置，包括主推进器1与辅助推进器2，通过操纵方向舵3改变航向进行自守卫航行，直至到达守卫内圆中，然后停机继续进行测量任务。

如附图3所示，为本发明的机动浮标在随机海浪干扰下的自守卫航行控制方法流程图，主要包括以下内容：

步骤1、根据机动浮标在随机海浪干扰下的纵摇、横摇、垂荡运动的时间序列，采用mar(multiautoregressivemodel)模型法计算机动浮标在上述三自由度上运动的交叉谱；

步骤2、根据机动浮标自身的模型参数通过matlab中的mss工具包计算出浮标在上述三自由度上运动的rao(responseamplitudeoperator)数据；

步骤3、根据步骤1中获得的交叉谱和步骤2中获得的rao数据计算得到机动浮标的海浪方向谱，三者之间的关系为：

式中，φij(ω)为通过机动浮标的时间序列获取的交叉谱；raoi(ω,θ)表示在第i个自由度上，波浪入射角为θ，频率为ω时的响应幅值算子，s(ω,θ)为海浪方向谱；

步骤4、根据机动浮标的实际航行情况，在给定的自守卫区域内，设置虚拟节点，每个虚拟节点中都包含此时刻的浪向角信息；

步骤5、根据机动浮标在随机海浪干扰下受到的波浪增阻与浪向角之间的函数关系，将机动浮标在自守卫航行时的能量损耗作为目标函数；

步骤6、采用a*算法对机动浮标在随机海浪干扰下自守卫航行时，从守卫区域边缘点到达给定守卫点进行路径规划，使得机动浮标在规划的航路上航行时能量损耗最低。

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，步骤1中mar模型建立方法为：

步骤1、假设{zk}是k维随机平稳的零均值时间序列，表示如下：

zk＝[z1z2...zk]^t

其中，zi(s·δt),s＝1,2,...,n，i＝1,2,...,k,其中δt为采样时间，令采样时间δt＝1，则zi(s·δt)＝zi(s)；

步骤2、由于所需要分析的是机动浮标三自由度的运动时间序列，在每个时间序列中，其前p个时间序列值对第s个时间序列进行建模，则有：

z(s)-a1z(s-1)-...-ap(s-p)＝ε(s)

其中，ai为k×k维的权重系数矩阵；ε(s)＝[ε1ε2...εk]^t表示均值为0，方差为σi的高斯白噪声序列，由于εi与z(s)的取值无关，因此ε(s)与ε(s-p)的协方差为零；由此可得白噪声向量的协方差矩阵如下：

式中，σij表示εi与εj之间的协方差；将上式改写如下：

令：v(s)＝[z(s-1)z(s-2)...z(s-p)]^t

λp＝[λ1λ2...λp](k×(k·p))

则mar模型的最终形式为：

z(s)＝λpv(s)+ε(s)

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，对上述权重系数矩阵求解时采用whittle递推法；由于mar模型的权值系数aj(j＝1,2,...,p)随着模型阶数p的变化而改变；因此，令：λpj＝(ap1ap1...apj)；j＝1,2,...,p，ri(i＝1,2,...,p)为相关矩阵，对于不同迟滞l相关矩阵中的第i,j个元素的计算公式表示为：

因此，递推公式为：

针对上述递推法求出的p项权重系数矩阵采用最优阶数判断准则来选取最有效的一组权重系数矩阵；该准则的描述为：对于长度为n的k维平稳随机序列{zk}，其某一步预报误差方差矩阵的估计量用dp来表示，最终的预报误差为其行列式的值：

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，步骤1中计算机动浮标运动交叉谱的方法包括以下步骤：

根据上述建立的mar模型，选择适当的阶数，就可以准确地使用前p个历史值准确预测所求摇荡时序的当前值。

步骤1、由交叉谱密度的基本定义得到白噪声的交叉谱密度

式中f为变换后的频率，由于当且仅当s＝t时，即等式右边只是将σst进行傅里叶变换，而σst是一个常数项，对于常数项的傅里叶变换仍然是其本身；可得

步骤2、考虑残差的影响，令aij(0)＝-δij，其中：

根据上式可得：

其中，εi(s)是σj(s)的线性变换，模型当维数退化到一维时，aij(p)为脉冲响应函数：

此时，原始信号经过脉冲响应函数的线性变化后，在其g频点的功率谱密度phh(g)为：

phh(g)＝|b(g)|²pyy(g),-∞＜g＜∞

式中，pyy(g)是y(t)的功率谱密度，为脉冲响应函数h(τ)的频率响应函数；

分析可知，多维与一维的情况相同，但是多维误差的交叉谱为复数，可知：

其中，α,β＝1,2,...k；p(f)即为所求的交叉谱；考虑到其矩阵形式表示为：

π＝a(f)p(f)a(f)^h

式中，π为由σst组成的矩阵；假定矩阵a(f)为非奇异矩阵，上式两端同时乘以a(f)的逆矩阵和a(f)^h(即a(f)的共轭转置矩阵)，可得：

p(f)＝a(f)^-1πa(f)^-h

综上，便可获得机动浮标三自由度摇荡运动的交叉谱，最终可得海浪方向谱。

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，机动浮标在随机海浪干扰下产生的波浪增阻f(β)与浪向即遭遇角β之间的函数关系为：

f(β)＝p0+p1β+...+pkβi^k

式中，(p0,...,pk)表示所需要求出的多项式系数，在一级海况、二级海况、三级海况下，采用最小二乘法对试验数据进行拟合仿真，求得的系数矩阵分别为：p＝[2.73×10^-10-1.22×10^-71.85×10^-4-7.44×10^-4]^t，p＝[-1.24×10^-98.85×10^-7-2.22×10^-40.027]^t，p＝[-1.13×10^-91.53×10^-6-4.92×10^-40.07]^t。

所述一种机动浮标在随机海浪干扰下的自守卫控制方法，步骤6具体包括以下步骤：

步骤6.1、设计代价函数，a*算法中的代价函数f(x)表示为：

f(x)＝g(x)+h(x)

式中，f(x)表示从起始节点n0开始到终止节点的最佳路径的代价总值；g(x)表示从起始节点n0开始到某一个节点x时的最佳路径代价值，h(x)表示从某一节点x到终止节点的实际最佳路径代价值，h(x)又被称为启发函数，当h(x)的值小于某一节点x到终止节点的实际距离时，可以得到最优解；

步骤6.2、根据波浪增阻与遭遇角之间的函数关系表达式，可得机动浮标航行至某一节点处的能量消耗值，设机动浮标的航速为v，起始节点的位置为(x0,y0)，在某一节点ni(xi,yi)节点处，航行至下一节点ni+1(xi+1,yi+1)的时间为t，则在这两个相邻节点的路径长度为l，则：

机动浮标在这段路径上航行的能量损失δploss(ni,ni+1)表示为：

把在整个规划的路径上损耗的总能量ploss作为代价函数，表示为：

步骤6.3、设计启发函数，选择欧几里得距离法作为机动浮标路径规划的启发方法，其对应的函数表达式表示为：

如附图4所示，为本发明用于机动浮标在随机海浪干扰下的自守卫航行控制方法的虚拟节点设置示意图，根据机动浮标在湖上和海上试验的实际情况，自守卫区域半径一般取50m左右，守卫内径在5m左右，而机动浮标自身的长度为1.35m。因此，本发明设置的虚拟节点是以自守卫内圆为中心点，在守卫区域内以5m的间距均匀设置虚拟节点。每个虚拟节点都包含此时的浪向角信息。其中，规定机动浮标开始自守卫航行时的点为起始节点n0，到达设置的守卫目标点为终止节点nn。

如附图5所示，为本发明用于机动浮标在随机海浪干扰下的自守卫航行控制方法的a*算法设计流程图，对其设计步骤详细描述如下：

本发明对于虚拟节点的设置与选取是依据机动浮标在随机海浪干扰下，产生的波浪增阻最小为准则；当机动浮标航行至某一节点时，其艏向与海浪方向之间的夹角即遭遇角的变化，影响波浪增阻的变化；通过航向控制器调整艏向来改变遭遇角的大小，进而使波浪增阻减少。此时选取多个虚拟节点，通过计算这些到达节点处的波浪增阻大小，通过比较来判断选取最优节点，进而完成能量最优路径规划。

步骤1、建立三个表，分别为open表、closed表、parent表，其中，open表用来存放当前根据海流和海浪方向信息规划的节点信息、closed表用来存放已经选取的节点信息、parent表用来存放父节点信息；

步骤2、对算法进行初始化设置，即将open表、closed表和parent表清空，然后设置起始节点n0、终止节点nn的位置信息；

步骤3、将起始节点n0放入open表中，令g(n0)＝0，则f＝h；

步骤4、判断open表是否为空，若open表为空，则所需规划的路径不存在，算法结束，若open表不为空，则执行步骤5；

步骤5、把open表中的代价函数值最小的节点定义为nbest，判断nbest是否为终止节点nn，若nbest＝nn，则表示当前路径规划完成，算法结束，如果nbest≠nn，则执行步骤6；

步骤6、计算open表中所存入的节点的代价函数值ploss，并比较相对应的代价函数值的大小，把代价函数值最小的节点即nbest从open表中赋给closed表中；

步骤7、判断是否满足能量损耗最低的约束条件，若满足在该节点的能量损耗最低，则计算该节点的代价值，并与实际代价值进行比较，将代价值最小的节点添加到open表和parent表中，返回步骤6，若不满足能量损耗最低的约束条件，则返回步骤4。