基于无迹卡尔曼滤波的超空泡航行体最优控制研究

2019-08-06张成举王聪曹伟王金强

兵工学报 2019年6期

张成举，王聪，曹伟，王金强

(哈尔滨工业大学航天学院, 黑龙江哈尔滨 150001)

0 引言

当水下航行体运动速度到达50 m/s时，航行体头部或整体会被空泡包裹，此时航行体称为超空泡航行体。超空泡航行体高速运动时会出现两种情况：一是航行体未被空泡全部包裹，航行体滑行力不为0；二是航行体被空泡全部包裹，航行体滑行力为0.海洋环境复杂多变，超空泡航行体受环境噪声和测量噪声的影响较大，因此超空泡航行体的运动控制具有极大的不确定性，所以超空泡航行体的运动控制是当今的研究热点。

在超空泡航行体运动控制方面，已有许多学者开展了研究，例如：Dzieiski等[1]建立了纵平面内4状态2自由度模型，对其基准问题开展了研究；Vanek等[2]考虑超空泡航行体运动过程中的空泡存在记忆效应，基于反馈线性化方法提出了双闭环控制系统，证明了该控制系统的跟踪控制有效性；Qiang等[3]改进了反步法，提出状态反馈控制，对超空泡航行体的运动稳定性开展了研究；Nguyen等[4]考虑空泡形态变化会对超空泡航行体运动有重要影响，提出反馈控制，对超空泡航行体运动稳定性开展了研究；Lü等[5]基于保代价理论提出自适应滑模控制器，针对超空泡航行体存在非匹配不确定性问题开展了控制研究；Ruzzence等[6]运用最优方程获取航行体的最优弹道，对航行体转弯运动和下潜运动开展了研究；Mao等[7]对超空泡航行体提出基于反步法的自适应控制器，设计尾舵效率计算方法，利用李雅普诺夫方法证明了该方法的有效性；Zhao等[8]基于超空泡航行体空泡突变特性，对航行体变深控制进行了研究；白涛等[9]提出基于混合卡尔曼滤波的变深运动控制算法，针对超空泡航行体的测量噪声干扰开展了研究，采用高斯白噪声测量了噪声，仿真结果表明该控制系统提高了系统的运动稳定性；庞爱平等[10]提出采用鲁棒控制算法来保证超空泡航行体各个状态变量在指定范围内变化，通过仿真证明该系统具有较好的稳定性；吕瑞等[11]针对超空泡航行体提出基于增益自适应变结构控制算法，经过仿真验证了该控制算法的有效性；王京华等[12]提出基于空泡记忆效应的超空泡航行体控制算法，通过仿真验证了该控制算法的有效性；李洋等[13]提出针对超空泡航行体的反演控制算法，通过仿真验证了该控制系统的稳定性。

根据以上文献可知，采用滤波器同时降低环境干扰与测量干扰的研究较少。本文针对超空泡航行体受到环境噪声和测量噪声的影响，设计基于无迹卡尔曼滤波器(UKF)的最优控制器，通过仿真分析验证了该控制器的有效性。

1 超空泡航行体动力学模型

以超空泡航行体纵平面内稳定问题为研究对象，通过对航行体建模受力进行分析，可知超空泡航行体主要受力如下：

1) 航行体所受重力G，重力方向指向地心。

2) 超空泡航行体发动机的推力T，推力作用方向与航行体纵轴重合。

3) 空化器转动过程中为航行体提供的升力[1]：

(1)

式中：ρ为水的密度；Rn为空化器半径；v为航行体航行速度；c0为空化器初始阻力系数；σ为航行体空化数；vy为航行体纵向速度；Lc为空化器顶端至航行体质心之间的距离；ωz为航行体俯仰角速度；δc为空化器转角。

4) 超空泡航行体航行过程中，其空泡区域分为全沾湿区、半沾湿区和空泡包裹区。航行体空泡区域划分如图1所示，其中δf为航行体尾翼转角。

图1 超空泡航行体空泡区域示意图Fig.1 Configuration of supercavitating vehicle’s cavitation

空化器转动会使空泡轴线偏移，重力作用也会导致空泡上漂。对空泡形态预测利用较多的是Logvinovich[14]提出的空泡形态半理论和半经验预测公式。本文采用此公式用于空泡形态预测，在文献[11]提出的空泡状态方程中，空泡截面为轴对称椭圆，假设空泡在任意位置均为圆形，t时刻空泡半径为

(2)

(3)

空泡受重力影响会导致空泡轴线发生偏移，在空化数σ<0.1且弗劳德数Fr较大的条件下，重力导致的空泡轴线偏移量为

(4)

式中：g为重力加速度；LRmax为空泡长度,LRmax=2Rn/(1.92/σ-3)；ζ=x/LRmax，x为空化器与空泡截面的距离。

空化器转动对空泡形态影响较大，并且航行体俯仰运动对空泡产生极大的非线性干扰，基于独立膨胀原理，空化器转动导致空泡轴线偏移量为

hc=0.82(1+σ)(θ+δc)Rn(0.46-σ+ζ)，

(5)

式中：θ为航行体俯仰角。

当超空泡航行体部分处于沾湿状态时，航行体尾部会与空泡之间发生接触和碰撞，航行体产生滑行力。根据细长体理论，将滑行力简化为工程中计算细长体浸入自由液面，采用如下考虑重力和空化器转动的滑行力计算公式[2]：

(6)

式中：

(7)

(8)

R为航行体半径，h′为沾湿深度，αp为航行体沾湿角，yc为空泡轴线与航行体轴线之间的纵向偏差，zc为空泡轴线与航行体轴线之间的侧向偏差，ω为航行体旋转角速度。

5) 尾舵包括俯仰舵和偏航舵，当尾翼处于沾湿状态时，尾舵才能发挥作用；航行体运动过程中，尾翼与空泡之间的关系处于不确定状态；假定航行体在纵平面内做小角度机动，此时可认为尾翼是一种楔形空化器，与空化器估算公式相同，则尾翼升力估算公式为

(9)

式中：η为尾翼升力与空化器升力的相似准数，这里取0.5；Lf为尾翼与航行体质心之间的距离。

建立动量定理和动量矩定理，可得超空泡航行体纵向运动模型矩阵：

(10)

(11)

I(t,τ)为尾舵沾湿率，

(12)

2 超空泡航行体最优控制器设计

采用精确线性化方法，利用非线性反馈变换实现输入输出的精确线性化[15]，首先将方程(10)式改写为

(13)

式中：

(14)

(15)

(16)

(17)

(18)

y1、y2为输出变量y的分量，n1、n2为充分光滑的标量函数。

对于此非线性多输入多输出(MIMO)系统，具有相对阶向量(λ1,λ2)=(2,2)，总相对阶[15]为λ1+λ2=4，这时不考虑内动态的稳定性，输出变量与输入的关系为

(19)

式中：Lχ、LMi(i=1,2)分别为向量χ(x)和M(x)的李导数。

通过反馈可得

(20)

将(20)式代入(19)式，可得

(21)

将(21)式精确线性化后，可得超空泡航行体状态方程：

(22)

(23)

式中：ε(t)为线性二次型最优控制器；k1、k2、k3、k4为控制参数，其选择应使两个解耦子系统对应特征多项式的所有根均位于左半复平面；h1、h2、h3、h4为向量h(t)的分量，则可得误差方程为

(24)

式中：

(25)

选取如下线性二次型性能指标离散泛函数：

(26)

式中：Q为对称正定状态加权矩阵；R为对称正定控制加权矩阵。

当方程(26)式取得最小值时，可得线性二次型最优控制反馈控制器为

ε(t)=-R-1ΘTPh(t)，

(27)

式中：P为满足Riccati方程(28)式的正定对称矩阵，

HTP+PH-PΘR-1ΘTP+Q=0.

(28)

综上所述，联立方程(20)式、方程(23)式和方程(27)式，可得基于精确线性化后的最优控制器表达式为

(29)

3 带有环境噪声和测量噪声的最优控制算法仿真分析

令：超空泡航行体初始纵向速度vy=-4 m/s，俯仰角θ=3°，航行深度y=-4 m,俯仰角速度ωz=7°/s；正弦跟踪方程为y(t)=-sint-5；航行体航行速度v=75 m/s，质量m=22 kg，重力加速度g=9.81 m/s2,空化器半径Rn=0.019 1 m，航行体柱段半径R=0.050 8 m，航行体总长L=1.8 m，空化数σ=0.03，水的密度ρ=998 kg/m3,Lc=17L/28,Lf=11L/28.

采用经输入输出精确线性化后的最优控制器对超空泡航行体进行位置跟踪控制；环境噪声和测量噪声幅值均为0.02 m/s的高斯白噪声[9]。考虑到较小的控制误差，经过多次仿真实验，为保证控制精度，选取控制参数k1=k2=k3=k4=20，

(30)

将Q和R值代入方程(29)式，满足P为正定矩阵。采用以上参数进行仿真，验证算法的有效性，仿真时间为20 s，仿真结果如图2～图8所示。

图2 航行体深度变化(无滤波器)Fig.2 Change of depth (without filter)

图3 航行体纵向速度变化(无滤波器)Fig.3 Change of vertical speed (without filter)

由图2～图8可知，超空泡航行体在环境噪声和测量噪声干扰下出现持续无规则振荡。由于环境噪声和测量噪声的干扰，航行体滑行力呈现无规则振荡，给航行体运动控制带来了困难。

图4 航行体俯仰角速度变化(无滤波器)Fig.4 Change of pitch angular velocity (without filter)

图5 航行体俯仰角变化(无滤波器)Fig.5 Change of pitch angle (without filter)

图6 航行体滑行力变化(无滤波器)Fig.6 Change of sliding force (without filter)

根据上述对超空泡航行体在环境噪声和测量噪声干扰下的运动状态分析可知，超空泡航行体的运动是极其不稳定的。因此，在超空泡运动控制中有必要加入卡尔曼滤波器，使航行体减少干扰噪声的影响。

图7 航行体空化器转角变化(无滤波器)Fig.7 Change of cavitor angle (without filter)

图8 航行体尾翼转角变化(无滤波器)Fig.8 Change of tail angle (without filter)

4 基于UKF的最优控制算法

本文针对以上超空泡航行体的运动不稳定性，提出基于UKF的最优控制算法，解决超空泡航行体的噪声干扰问题。UKF采用卡尔曼线性滤波框架，对于一步预测方程，使用如下无迹变换处理均值和协方差的非线性传递问题：

(31)

式中：x(k+1)为k+1时刻的系统状态；f(x(k),W(k))为非线性状态方程函数，W(k)为环境噪声；y(k)为观测变量；h(x(k),V(k))为非线性观测方程函数，V(k)为测量噪声；W(t)为协方差S1的连续时间高斯白噪声；V(t)为协方差S2的连续时间高斯白噪声。

具体处理步骤如下：

1) 获取超空泡航行体的一组采样点，设置如下初值：

x(i)(k|k)=
[(k|k)(k|k)+T(k|k)-T]，

(32)

式中：x(i)(k|k)为k时刻超空泡航行体状态值，i为第i个采样点，i=1,2,…,2n，n为状态维数；(k|k)为k时刻状态估计值；为半正定矩阵，D(k|k)为状态x(k|k)的方差，λ为缩放比例参数，λ=α2(n+a)-n，用于减小总的预测误差；α为控制采样点分布状态的参数，用于调整Sigma点与均值点的距离。

2) 计算2n个Sigma点集的一步预测：

x(i)(k+1|k)=f[k,x(i)(k|k)]，

(33)

式中：x(i)(k+1|k)为超空泡航行体k+1时刻的预测状态；f[k,x(i)(k|k)]为k时刻的非线性状态方程。

3) 预测超空泡航行体的系统状态量与方差：

(34)

式中：wm为采样点均值的权值；wc为采样点协方差的权值。

采样点均值的权值和协方差的权值计算公式如下：

(35)

4) 根据一步预测值，再次运用无迹变换产生新的Sigma点集：

(36)

5) 量测更新方程，得到Sigma点集的观测预测值，通过加权求和得到超空泡航行体系统预测的均值、状态量和观测量之间的协方差、新息序列的协方差分别为

(37)

6) 计算卡尔曼增益矩阵：

(38)

7) 计算系统的状态更新和协方差更新：

k(k+1|k+1)=k(k+1|k)+Kk(yk-k)，
D(k+1|k+1)=

(39)

5 基于UKF的最优控制仿真分析

令：超空泡航行体初始纵向速度vy=-4 m/s，俯仰角θ=3°，航行深度y=-4 m,俯仰角速度ωz=7°/s；α通常选取一个较小的正实数(一般10-4≤α≤1)，α=0.01；a为待选参数，需要保证矩阵(n+λ)D(k|k)为半正定矩阵；κ=0；β用于合并方程中高阶项的动差，调节β可改善方差的精度，对于高斯分布，β=2时最优。选取的控制器参数k1=k2=k3=k4=20，Q和R与第3节相同，对选取的深度信号y(t)=-sint-5进行深度跟踪，仿真时间20 s时得到仿真结果如图9～图15所示。

图9 航行体深度变化(含滤波器)Fig.9 Change of depth (filter)

由图9～图14可知，通过运用UKF对噪声干扰处理，超空泡航行体处于稳定运动状态。

由图9可知，经过UKF处理后的跟踪信号与理想信号相差较小，误差能够控制在0.1 m内，而且未出现跟踪延迟现象。

由图10可知，超空泡航行体俯仰角速度急剧反向转变，在4 s时趋于稳定状态，经过滤波器处理后俯仰角速度跟踪误差控制较小。

图10 航行体俯仰角速度变化(含滤波器)Fig.10 Change of pitch angular velocity (filter)

图11 航行体纵向速度变化(含滤波器)Fig.11 Change of vertical speed (filter)

由图11可知，超空泡航行体纵向速度变化稳定，跟踪误差在0.1 m/s范围内。

由图12可知，超空泡航行体俯仰角逐渐较小，变化率逐渐减小，在4 s时航行体俯仰角趋于稳定，经过UKF处理后的俯仰角跟踪误差控制较小。

图12 航行体俯仰角变化(含滤波器)Fig.12 Change of pitch angle (filter)

图13 航行体滑行力变化(含滤波器)Fig.13 Change of sliding force (filter)

图14 航行体空化器转角变化(含滤波器)Fig.14 Change of cavitor angle (filter)

由图13可知，超空泡航行体起始状态时航行体包裹面积较大，滑行力较小。开始运动后，滑行力呈现先增大、后减小现象，之后航行体滑行力减小为0，此时航行体除空化器外，均为空泡包裹，经过UKF处理后的滑行力信号与理想状态相差无几。在航行体滑行力稳定之前，误差控制在0.1 kN范围内，但是从起始状态到稳定状态，经过UKF处理后的滑行力信号仍然比理想信号延迟。

由图14和图15可知，航行体空化器转角与尾翼转角偏转范围较小，与理想信号较为接近。