基于CS-Jounce模型的船摇滤波与预报方法

戴正旭,杜昌平,郑 耀,陈嘉鸿

(1.浙江大学航空航天学院,浙江 杭州 310027;2.中国卫星海上测控部,江苏 江阴 214431)

0 引 言

海上动平台测量需要获取船舶姿态数据,进行姿态修正后才能得到准确的目标位置信息[1],对实测船摇数据滤波去噪能提高测量数据精度。跟踪目标时,为了提高跟踪精度,克服高动态工况下动态滞后误差较大的问题,需要计算前馈值,输入天线伺服跟踪环路[2],此时需要计算船摇数据角度量、角速度量及角加速度量的预报值。故研究船摇数据特性及其滤波和预报方法,对提高海上综合测量精度有较强的意义。

根据工程应用需求,船摇数据实时滤波及预报需满足以下要求:①数据积累过程短,达到满足工作条件的积累点数尽量少;②计算量较小,能满足实时计算的要求;③能同时适用于横摇、纵摇、艏摇的数据滤波及预报;④能同时给出船摇角度数据、角速度及角加速度数据的滤波及预测值;⑤预报精度尽可能高;⑥鲁棒性较好,不易受野值影响。

国内外针对船舶姿态预报的研究成果很多,长自回归(autoregressive,AR)模型[3]在确定阶数的条件下精度满足要求,但有积累点数较多、不同海况或船舶工况下需调整参数的问题;卡尔曼滤波在船摇预报中的研究目前大多基于船舶水动力学模型[4-9],缺点是建模复杂,风浪的影响无法准确建模,且一旦动力学参数发生变化时,也会导致滤波精度下降;基于神经网络的方法[10-13]具有较高的预报精度,但只针对横摇数据进行研究,对艏摇和纵摇数据适用情况未有涉及。基于信号分解和AR模型综合预报的方法[14]提高了船摇预报时长,但实时性能不能满足要求。此外以上所有方法均未涉及船摇数据的滤波问题,以及船摇角速度和角加速度的估计和预报问题。

文中将基于船摇数据特性分析,引入机动目标跟踪理论,提出基于最大相关熵的自适应“当前”统计Jounce模型(“current” statistical Jounce model,CS-Jounce),改进及创新点包括:①针对目前Jerk模型阶数不够的问题,推导了更高一阶的Jounce模型;②结合自适应方法和最大相关熵理论,提高了模型的跟踪稳定性和鲁棒性;③无需对数据进行特殊预处理,同时适用于船舶的横摇、纵摇及航向数据预报;④可同时解决船摇数据的滤波、预报及角速度、角加速度的估计和预报问题。

1 船摇数据特性分析

文中基于某测量船姿态数据进行分析,3级海况,采样频率20 Hz,数据点数8 000点,时长400 s,包括横摇、纵摇和航向。

1.1 时域分析

对信号的传统分析方法是波形分析,信号的特性首先表现为它的时间特性。从原始数据图(见图1)可以看出,船摇数据具有明显的非平稳性和周期性,横摇和纵摇数据均值在零附近,航向数据均值可能较大,且由于船舶航行中可能改变航向,除有局部短周期外还可能有长趋势非周期项。

图1 原始船摇数据Fig.1 Raw ship-swaying data

1.2 频域分析

图2给出了船摇数据的频谱特性分析图,因航向数据均值非零,低频部分均出现了峰值,掩盖了其他频谱特征部分,对此文中采用一次差分的方法消除低频特征后再研究(见图2(d))。从图2中可以看出,船摇数据均存在明显的频率峰值,但非单峰分布,横摇、纵摇、航向数据频谱分布带相似,集中在小于0.6 Hz的低频波段。

图2 船摇数据频谱特性图Fig.2 Spectrum of ship-swaying data

1.3 噪声特性分析

常用的测量数据噪声分析方法有时域、频域及时频分析方法。

时域分析中常用的有变量差分法和最小二乘拟合残差法[15]。但两者都有局限性:变量差分法对于随机误差序列非白噪声的数据统计结果会失真;最小二乘拟合残差法在未加权情况下统计结果容易受到野值的影响,且拟合使用的数据点数及多项式阶数确定比较繁琐。

基于频域的滤波分析法,是基于傅里叶变换,将有用信号和噪声信号在频域进行分离,从而达到去噪的效果,缺陷在于,实际情况噪声频谱几乎都是分布在整个频域内,进行噪声平滑的同时,也必定平滑了非平稳信号的突变点。时频分析中的小波方法能克服频域滤波的缺点,但实际应用中需综合权衡分解阶数、阈值及小波核的选取。考虑到船摇数据基于船舶自身的惯性,不存在很多突变点,文中采用低通滤波的方法分离并分析船摇数据噪声。

设计巴特沃斯低通滤波器,通带截止频率0.6 Hz,阻带截止频率4 Hz,表1给出了各数据噪声均方差统计结果,图3分析了随机误差特性,其频谱图非均匀分布,具有隐周期成分,故船摇数据环境非高斯白噪声。

表1 船摇数据噪声均方差统计结果Table 1 Statistical results of the mean variance of ship-swaying data noise

图3 噪声频谱特性图Fig.3 Noise spectrum diagram

1.4 船摇数据高阶量分析

去除噪声后船摇数据高阶量信息可以通过差分的方法获得。图4横摇为例,绘制了滤波后数据的一阶至四阶差分曲线图。

图4 横摇高阶数据示意图Fig.4 High order data of rolling data

从图4中可以看出,船摇数据四阶差分前一直具有明显的周期性,对这样的数据使用Jerk模型明显阶数不足。四次差分后加加加速度量(jounce)周期性显著减弱,故船摇过程可以看成是jounce噪声量驱动过程。同时需注意到个别段落存在jounce量振荡加剧的现象。

表2给出了船摇数据Jounce量的统计结果。

表2 船摇数据jounce统计结果Table 2 Statistical results of ship-swaying data jounce

2 CS-Jounce模型

第1节分析指出船摇数据具有以下特征:

(1) 数据环境非高斯白噪声;

(2) 可以看成是jounce噪声量驱动过程,同时必须考虑个别时间段jounce振荡加剧的现象。

故做短时间预报时,可将船摇数据看作连续机动的过程,避免复杂的船舶姿态动力学建模。

常用的目标机动模型有[16]:匀速直线运动模型,匀加速直线运动模型,一阶时间相关模型,增加阶数的时间相关模型(Jerk)[17],“当前”统计模型等。在此基础上有文献[18]借鉴“当前”统计思想,提出了CS-Jerk模型。文献[19]分析了Jerk模型的缺陷,指出其跟踪高动态目标时稳态误差不为零,证明CS-Jerk可克服此问题。

从船摇数据特性可以看出使用Jerk算法仍有阶数不足的问题,故文中基于CS-Jerk模型思路,增加一阶,提出CS-Jounce模型。

2.1 CS-Jounce模型

船摇数据jounce量可认为服从非零均值的时间相关过程,将jounce上一时刻的一步预报值作为当前时刻jounce的均值,运动模型为

目标的连续系统状态方程为

将连续方程离散化后得目标的离散系统状态方程为

量测方程为

Z(k)=H(k)X(k)+V(k)

当采样周期为T时,状态转移矩阵为

其中

T=tk+1|-tk

o1=(6e-αT-6+6αT-3α2T2+α3T3)/(6α4)

p1=(2-2αT+α2T2-2e-αT)/(2α3)

q1=(e-αT-1+αT)/α2

r1=(1-e-αT)/α

s1=-e-αT

当αT足够小时,矩阵退化为

其中

q11=(18+36αT-36α2T2+24α3T3-12α4T4+

21α5T5/5-α6T6+α7T7/7-72αTe-αT-

12α3T3e-αT-18e-2αT)/36α9

q12=q21=(6-12αT+12α2T2-8α3T3+7α4T4/2-

α5T5+α6T6/6-12e-αT+12αTe-αT-

6α2T2e-αT+2α3T3e-αT+6e-2αT)/12a8

q13=q31=(3+6αT-6α2T2+3α3T3-α4T4+

α5T5/5-12αTe-αT-α3T3e-αT-3e-2αT)/6α7

q14=q41=(3-6αT+3α2T2-α3T3+α4T4/4-

6e-αT+6αTe-αT+α3T3e-αT+3e-2αT)/6α6

q15=q51=(3-6αTe-αT-α3T3e-αT-3e-2αT)/6α5

q22=(-3+2αT-2α2T2+4α3T3/3-α4T4/2+

α5T5/10+4e-αT+2α2T2e-αT-e-2αT)/2α7

q23=q32=(1-2αT+2α2T2-α3T3+α4T4/4-

2e-αT+2αTe-αT-α2T2e-αT+e-2αT)/2α6

q24=q42=(-3+2αT-α2T2+α3T3/3+4e-αT+

α2T2e-αT-e-2αT)/2α5

q25=q52=(1-2e-αT-α2T2e-αT+e-2αT)/2α4

q33=(1+2αT-2α2T2+2α3T3/3-4αTe-αT-e-2αT)/2a5

q34=q43=(1-2αT+α2T2-2e-αT+

2αTe-αT+e-2αT)/2α4

q35=q53=(1-2αTe-αT-e-2αT)/(2α3)

q44=(-3+2αT+4e-αT-e-2αT)/2α3

q45=q54=(1-2e-αT+e-2αT)/2α2

q55=(1-e-2αT)/2α

当αT足够小时,Q(k)退化为

Q(k)=

jounce的预测方差为

假设umax为最大机动jounce量,则jounce方差为

2.2 jounce自适应模型设计

CS模型缺点在于对目标机动频率、目标机动最大加速度两个参数依赖性较高。针对目标机动最大加速度自适应问题有很多改进方法,函数模糊隶属函数[20-23]是自适应调整中采用比较多。但相比于阶跃特征的信号,船舶姿态运动jounce量是在最大值范围内绕零值上下波动的,按CS自适应理论,处于零值附近时目标机动范围可能性更广,对应方差应该最大,但若按上述论文思路,将每个时刻的Jounce值都调整到[(4-π)umax/4,umax]或[u-max,(4-π)u-max/4]范围内显然是不可行的。

文中需要解决的问题是当船摇jounce量出现短时剧烈震荡时,其范围可能超出事前统计最大值umaxt,此时需要动态调整umax。以下借用高斯函数模糊隶属函数思想设计函数,原则为,设计最大安全机动范围uth>umaxt,正常情况下计算值umax≈umaxt,当u>umaxt时,动态计算保证umax>u的公式为

M函数借助了以uth为中心,δ为方差的正态分布函数,μ值保证1-μ≤M≤1,从而控制umax计算最小值。

根据表2中横摇jounce量的统计特性可知:umaxt=0.818,设计参数,δ=1.5,uth=4,μ=0.8。从隶属函数图5中可以看出,-1≤u≤1时umax变化平缓,|u|≥1时,umax能迅速调整适应更剧烈的机动状态。

图5 隶属函数示意图Fig.5 Membership function

3 最大相关熵卡尔曼滤波器

实际的船摇数据在数据的收集、传输或处理过程中经常受到一些随机误差的影响,产生各种噪音,它们极少服从正态分布,有的幅值较大,还有的是奇异点。

卡尔曼滤波器为基于最小方差理论针对线性系统的状态估计方法,因其实现简单且在高斯环境下获得最优解,在现实中被大量应用。但当系统环境非高斯白噪声时,滤波效果会变差[24],尤其数据含脉冲性质的噪声影响时,系统的鲁棒性欠佳。

针对非高斯白噪声系统的状态估计,常用方法有状态扩展、量测差分[25]。但对于重尾分布的噪声,很难基于一步相关模型去描述,同时也需要比较大的计算量,很难在实时应用。多模型方法[26]也是针对非高斯白噪声系统的手段,但主要缺点是多个模型并行运行,计算量较大。

最大相关熵准则(maximum correntropy criterion,MCC)作为非高斯信号处理的最佳鲁棒准则,近年来已成功应用于自适应滤波领域[27-28]。文献[29]结合最大相关熵准则和不动点迭代算法,提出了最大相关熵卡尔曼滤波器(maximum correntropy Kalman filter,MCKF)。理论和实践均表明,该滤波器对非高斯噪声具有良好的鲁棒性。

文中结合“当前”概率模型,总结算法流程概况如下:

步骤2根据状态预测和方差预测获得一步预报值，即

P(k|k-1)=

F(T)P(k-1|k-1)FT(k-1)+Q(k-1)

计算P(k|k-1)的Cholesky分解获得Bp(k|k-1);计算观测噪声阵R的Cholesky分解获得Br(k)。

(1)

其中

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(12)

步骤5比对当前状态估计量和步骤4迭代状态量,若满足(13)式的迭代收敛条件则终止迭代,进行步骤6,否则t+1→t,继续步骤4迭代计算。

(13)

步骤6更新最终方差阵,k+1→k,回到步骤2 进行下一点计算。

4 结果分析

4.1 参数设置

文中方法使用参数设置及原则如下:

(1) 目标机动常数的倒数α。根据表2中jounce频率分布范围,α取值范围为0.5～10,结果差别不大,根据频谱特征最集中段,分析时使用α=3。

(2) 观测噪声R。根据表1中噪声统计结果可设对应的观测噪声量。

(3) jounce机动最大值自适应参数。详见第2.2节。

(4) 相关熵参数设置。包括核函数带宽σ和迭代收敛条件ε。文献[29]中指出,对于非脉冲噪声,σ选取决定其性能是否能超过普通卡尔曼滤波器。以下分析时,综合考虑效果和性能,选择σ=3,ε=10-6。

4.2 实际效果

衡量滤波结果好坏的指标通常有均方根误差和平滑度,其中平滑度可以限制过拟合,但对于船摇数据平滑度指标设置比较困难,预测误差也可以反映模型准确度,限制过拟合,故文中结合滤波均方根误差和三步预测误差值的均方根误差综合考量。

图6以横摇数据为例给出了滤波效果,表3给出了滤波和预报误差统计结果。

图6 横摇数据滤波及预报效果图Fig.6 Filtering and forecasting effect of rolling data

表3 滤波及三步预报误差均方根Table 3 Mean square error of filter and three-step prediction

可以得出以下结论:

(1) 该算法能很好地跟踪船摇数据,滤波误差中不存在明显趋势量;滤波后数据平滑度较高,能有效滤除船摇数据噪声。

(2) 滤波误差统计结果与表1分析结果一致性较好。

(3) 表3的三步预报误差统计结果表明,该算法能实现船摇数据的准确预报,均方根误差小于5″。

(4) 表3的统计结果表明,该方法能准确预报船摇角速度及角加速度,角速度预报均方根误差约为0.006(°)/s,角加速度预报均方根误差小于0.4(°)/s2。

(5) 图6(d)表明,该方法鲁棒性较好,尤其针对脉冲噪声数据,不会出现普通卡尔曼滤波器的精度变差的现象。

5 结 论

本文通过分析航天测量船的船摇数据特性,引入机动目标跟踪理论进行船摇滤波及预报,针对现有模型阶数不够高的问题导出了四阶的CS-Jounce模型,并给出一种机动参数自适应模型。针对船摇数据环境非高斯白噪声的问题,引入了最大相关熵卡尔曼滤波器。实测数据分析证明,该方法能很好地跟踪船摇数据,滤波效果较理想,预报精度较高,且能很好地估计并预报船摇角速度及角加速度。