戴正旭,杜昌平,鄭 耀,陳嘉鴻

(1.浙江大學(xué)航空航天學(xué)院,浙江 杭州 310027;2.中國(guó)衛(wèi)星海上測(cè)控部,江蘇 江陰 214431)

0 引 言

海上動(dòng)平臺(tái)測(cè)量需要獲取船舶姿態(tài)數(shù)據(jù),進(jìn)行姿態(tài)修正后才能得到準(zhǔn)確的目標(biāo)位置信息[1],對(duì)實(shí)測(cè)船搖數(shù)據(jù)濾波去噪能提高測(cè)量數(shù)據(jù)精度。跟蹤目標(biāo)時(shí),為了提高跟蹤精度,克服高動(dòng)態(tài)工況下動(dòng)態(tài)滯后誤差較大的問(wèn)題,需要計(jì)算前饋值,輸入天線伺服跟蹤環(huán)路[2],此時(shí)需要計(jì)算船搖數(shù)據(jù)角度量、角速度量及角加速度量的預(yù)報(bào)值。故研究船搖數(shù)據(jù)特性及其濾波和預(yù)報(bào)方法,對(duì)提高海上綜合測(cè)量精度有較強(qiáng)的意義。

根據(jù)工程應(yīng)用需求,船搖數(shù)據(jù)實(shí)時(shí)濾波及預(yù)報(bào)需滿(mǎn)足以下要求:①數(shù)據(jù)積累過(guò)程短,達(dá)到滿(mǎn)足工作條件的積累點(diǎn)數(shù)盡量少;②計(jì)算量較小,能滿(mǎn)足實(shí)時(shí)計(jì)算的要求;③能同時(shí)適用于橫搖、縱搖、艏搖的數(shù)據(jù)濾波及預(yù)報(bào);④能同時(shí)給出船搖角度數(shù)據(jù)、角速度及角加速度數(shù)據(jù)的濾波及預(yù)測(cè)值;⑤預(yù)報(bào)精度盡可能高;⑥魯棒性較好,不易受野值影響。

國(guó)內(nèi)外針對(duì)船舶姿態(tài)預(yù)報(bào)的研究成果很多,長(zhǎng)自回歸(autoregressive,AR)模型[3]在確定階數(shù)的條件下精度滿(mǎn)足要求,但有積累點(diǎn)數(shù)較多、不同海況或船舶工況下需調(diào)整參數(shù)的問(wèn)題;卡爾曼濾波在船搖預(yù)報(bào)中的研究目前大多基于船舶水動(dòng)力學(xué)模型[4-9],缺點(diǎn)是建模復(fù)雜,風(fēng)浪的影響無(wú)法準(zhǔn)確建模,且一旦動(dòng)力學(xué)參數(shù)發(fā)生變化時(shí),也會(huì)導(dǎo)致濾波精度下降;基于神經(jīng)網(wǎng)絡(luò)的方法[10-13]具有較高的預(yù)報(bào)精度,但只針對(duì)橫搖數(shù)據(jù)進(jìn)行研究,對(duì)艏搖和縱搖數(shù)據(jù)適用情況未有涉及?；谛盘?hào)分解和AR模型綜合預(yù)報(bào)的方法[14]提高了船搖預(yù)報(bào)時(shí)長(zhǎng),但實(shí)時(shí)性能不能滿(mǎn)足要求。此外以上所有方法均未涉及船搖數(shù)據(jù)的濾波問(wèn)題,以及船搖角速度和角加速度的估計(jì)和預(yù)報(bào)問(wèn)題。

文中將基于船搖數(shù)據(jù)特性分析,引入機(jī)動(dòng)目標(biāo)跟蹤理論,提出基于最大相關(guān)熵的自適應(yīng)“當(dāng)前”統(tǒng)計(jì)Jounce模型(“current” statistical Jounce model,CS-Jounce),改進(jìn)及創(chuàng)新點(diǎn)包括:①針對(duì)目前Jerk模型階數(shù)不夠的問(wèn)題,推導(dǎo)了更高一階的Jounce模型;②結(jié)合自適應(yīng)方法和最大相關(guān)熵理論,提高了模型的跟蹤穩(wěn)定性和魯棒性;③無(wú)需對(duì)數(shù)據(jù)進(jìn)行特殊預(yù)處理,同時(shí)適用于船舶的橫搖、縱搖及航向數(shù)據(jù)預(yù)報(bào);④可同時(shí)解決船搖數(shù)據(jù)的濾波、預(yù)報(bào)及角速度、角加速度的估計(jì)和預(yù)報(bào)問(wèn)題。

1 船搖數(shù)據(jù)特性分析

文中基于某測(cè)量船姿態(tài)數(shù)據(jù)進(jìn)行分析,3級(jí)海況,采樣頻率20 Hz,數(shù)據(jù)點(diǎn)數(shù)8 000點(diǎn),時(shí)長(zhǎng)400 s,包括橫搖、縱搖和航向。

1.1 時(shí)域分析

對(duì)信號(hào)的傳統(tǒng)分析方法是波形分析,信號(hào)的特性首先表現(xiàn)為它的時(shí)間特性。從原始數(shù)據(jù)圖(見(jiàn)圖1)可以看出,船搖數(shù)據(jù)具有明顯的非平穩(wěn)性和周期性,橫搖和縱搖數(shù)據(jù)均值在零附近,航向數(shù)據(jù)均值可能較大,且由于船舶航行中可能改變航向,除有局部短周期外還可能有長(zhǎng)趨勢(shì)非周期項(xiàng)。

圖1 原始船搖數(shù)據(jù)Fig.1 Raw ship-swaying data

1.2 頻域分析

圖2給出了船搖數(shù)據(jù)的頻譜特性分析圖,因航向數(shù)據(jù)均值非零,低頻部分均出現(xiàn)了峰值,掩蓋了其他頻譜特征部分,對(duì)此文中采用一次差分的方法消除低頻特征后再研究(見(jiàn)圖2(d))。從圖2中可以看出,船搖數(shù)據(jù)均存在明顯的頻率峰值,但非單峰分布,橫搖、縱搖、航向數(shù)據(jù)頻譜分布帶相似,集中在小于0.6 Hz的低頻波段。

圖2 船搖數(shù)據(jù)頻譜特性圖Fig.2 Spectrum of ship-swaying data

1.3 噪聲特性分析

常用的測(cè)量數(shù)據(jù)噪聲分析方法有時(shí)域、頻域及時(shí)頻分析方法。

時(shí)域分析中常用的有變量差分法和最小二乘擬合殘差法[15]。但兩者都有局限性:變量差分法對(duì)于隨機(jī)誤差序列非白噪聲的數(shù)據(jù)統(tǒng)計(jì)結(jié)果會(huì)失真;最小二乘擬合殘差法在未加權(quán)情況下統(tǒng)計(jì)結(jié)果容易受到野值的影響,且擬合使用的數(shù)據(jù)點(diǎn)數(shù)及多項(xiàng)式階數(shù)確定比較繁瑣。

基于頻域的濾波分析法,是基于傅里葉變換,將有用信號(hào)和噪聲信號(hào)在頻域進(jìn)行分離,從而達(dá)到去噪的效果,缺陷在于,實(shí)際情況噪聲頻譜幾乎都是分布在整個(gè)頻域內(nèi),進(jìn)行噪聲平滑的同時(shí),也必定平滑了非平穩(wěn)信號(hào)的突變點(diǎn)。時(shí)頻分析中的小波方法能克服頻域?yàn)V波的缺點(diǎn),但實(shí)際應(yīng)用中需綜合權(quán)衡分解階數(shù)、閾值及小波核的選取?？紤]到船搖數(shù)據(jù)基于船舶自身的慣性,不存在很多突變點(diǎn),文中采用低通濾波的方法分離并分析船搖數(shù)據(jù)噪聲。

設(shè)計(jì)巴特沃斯低通濾波器,通帶截止頻率0.6 Hz,阻帶截止頻率4 Hz,表1給出了各數(shù)據(jù)噪聲均方差統(tǒng)計(jì)結(jié)果,圖3分析了隨機(jī)誤差特性,其頻譜圖非均勻分布,具有隱周期成分,故船搖數(shù)據(jù)環(huán)境非高斯白噪聲。

表1 船搖數(shù)據(jù)噪聲均方差統(tǒng)計(jì)結(jié)果Table 1 Statistical results of the mean variance of ship-swaying data noise

圖3 噪聲頻譜特性圖Fig.3 Noise spectrum diagram

1.4 船搖數(shù)據(jù)高階量分析

去除噪聲后船搖數(shù)據(jù)高階量信息可以通過(guò)差分的方法獲得。圖4橫搖為例,繪制了濾波后數(shù)據(jù)的一階至四階差分曲線圖。

圖4 橫搖高階數(shù)據(jù)示意圖Fig.4 High order data of rolling data

從圖4中可以看出,船搖數(shù)據(jù)四階差分前一直具有明顯的周期性,對(duì)這樣的數(shù)據(jù)使用Jerk模型明顯階數(shù)不足。四次差分后加加加速度量(jounce)周期性顯著減弱,故船搖過(guò)程可以看成是jounce噪聲量驅(qū)動(dòng)過(guò)程。同時(shí)需注意到個(gè)別段落存在jounce量振蕩加劇的現(xiàn)象。

表2給出了船搖數(shù)據(jù)Jounce量的統(tǒng)計(jì)結(jié)果。

表2 船搖數(shù)據(jù)jounce統(tǒng)計(jì)結(jié)果Table 2 Statistical results of ship-swaying data jounce

2 CS-Jounce模型

第1節(jié)分析指出船搖數(shù)據(jù)具有以下特征:

(1) 數(shù)據(jù)環(huán)境非高斯白噪聲;

(2) 可以看成是jounce噪聲量驅(qū)動(dòng)過(guò)程,同時(shí)必須考慮個(gè)別時(shí)間段jounce振蕩加劇的現(xiàn)象。

故做短時(shí)間預(yù)報(bào)時(shí),可將船搖數(shù)據(jù)看作連續(xù)機(jī)動(dòng)的過(guò)程,避免復(fù)雜的船舶姿態(tài)動(dòng)力學(xué)建模。

常用的目標(biāo)機(jī)動(dòng)模型有[16]:勻速直線運(yùn)動(dòng)模型,勻加速直線運(yùn)動(dòng)模型,一階時(shí)間相關(guān)模型,增加階數(shù)的時(shí)間相關(guān)模型(Jerk)[17],“當(dāng)前”統(tǒng)計(jì)模型等。在此基礎(chǔ)上有文獻(xiàn)[18]借鑒“當(dāng)前”統(tǒng)計(jì)思想,提出了CS-Jerk模型。文獻(xiàn)[19]分析了Jerk模型的缺陷,指出其跟蹤高動(dòng)態(tài)目標(biāo)時(shí)穩(wěn)態(tài)誤差不為零,證明CS-Jerk可克服此問(wèn)題。

從船搖數(shù)據(jù)特性可以看出使用Jerk算法仍有階數(shù)不足的問(wèn)題,故文中基于CS-Jerk模型思路,增加一階,提出CS-Jounce模型。

2.1 CS-Jounce模型

船搖數(shù)據(jù)jounce量可認(rèn)為服從非零均值的時(shí)間相關(guān)過(guò)程,將jounce上一時(shí)刻的一步預(yù)報(bào)值作為當(dāng)前時(shí)刻jounce的均值,運(yùn)動(dòng)模型為

目標(biāo)的連續(xù)系統(tǒng)狀態(tài)方程為

將連續(xù)方程離散化后得目標(biāo)的離散系統(tǒng)狀態(tài)方程為

量測(cè)方程為

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

當(dāng)采樣周期為T(mén)時(shí),狀態(tài)轉(zhuǎn)移矩陣為

其中

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

當(dāng)αT足夠小時(shí),矩陣退化為

其中

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α

當(dāng)αT足夠小時(shí),Q(k)退化為

Q(k)=

jounce的預(yù)測(cè)方差為

假設(shè)umax為最大機(jī)動(dòng)jounce量,則jounce方差為

2.2 jounce自適應(yīng)模型設(shè)計(jì)

CS模型缺點(diǎn)在于對(duì)目標(biāo)機(jī)動(dòng)頻率、目標(biāo)機(jī)動(dòng)最大加速度兩個(gè)參數(shù)依賴(lài)性較高。針對(duì)目標(biāo)機(jī)動(dòng)最大加速度自適應(yīng)問(wèn)題有很多改進(jìn)方法,函數(shù)模糊隸屬函數(shù)[20-23]是自適應(yīng)調(diào)整中采用比較多。但相比于階躍特征的信號(hào),船舶姿態(tài)運(yùn)動(dòng)jounce量是在最大值范圍內(nèi)繞零值上下波動(dòng)的,按CS自適應(yīng)理論,處于零值附近時(shí)目標(biāo)機(jī)動(dòng)范圍可能性更廣,對(duì)應(yīng)方差應(yīng)該最大,但若按上述論文思路,將每個(gè)時(shí)刻的Jounce值都調(diào)整到[(4-π)umax/4,umax]或[u-max,(4-π)u-max/4]范圍內(nèi)顯然是不可行的。

文中需要解決的問(wèn)題是當(dāng)船搖jounce量出現(xiàn)短時(shí)劇烈震蕩時(shí),其范圍可能超出事前統(tǒng)計(jì)最大值umaxt,此時(shí)需要?jiǎng)討B(tài)調(diào)整umax。以下借用高斯函數(shù)模糊隸屬函數(shù)思想設(shè)計(jì)函數(shù),原則為,設(shè)計(jì)最大安全機(jī)動(dòng)范圍uth>umaxt,正常情況下計(jì)算值umax≈umaxt,當(dāng)u>umaxt時(shí),動(dòng)態(tài)計(jì)算保證umax>u的公式為

M函數(shù)借助了以u(píng)th為中心,δ為方差的正態(tài)分布函數(shù),μ值保證1-μ≤M≤1,從而控制umax計(jì)算最小值。

根據(jù)表2中橫搖jounce量的統(tǒng)計(jì)特性可知:umaxt=0.818,設(shè)計(jì)參數(shù),δ=1.5,uth=4,μ=0.8。從隸屬函數(shù)圖5中可以看出,-1≤u≤1時(shí)umax變化平緩,|u|≥1時(shí),umax能迅速調(diào)整適應(yīng)更劇烈的機(jī)動(dòng)狀態(tài)。

圖5 隸屬函數(shù)示意圖Fig.5 Membership function

3 最大相關(guān)熵卡爾曼濾波器

實(shí)際的船搖數(shù)據(jù)在數(shù)據(jù)的收集、傳輸或處理過(guò)程中經(jīng)常受到一些隨機(jī)誤差的影響,產(chǎn)生各種噪音,它們極少服從正態(tài)分布,有的幅值較大,還有的是奇異點(diǎn)。

卡爾曼濾波器為基于最小方差理論針對(duì)線性系統(tǒng)的狀態(tài)估計(jì)方法,因其實(shí)現(xiàn)簡(jiǎn)單且在高斯環(huán)境下獲得最優(yōu)解,在現(xiàn)實(shí)中被大量應(yīng)用。但當(dāng)系統(tǒng)環(huán)境非高斯白噪聲時(shí),濾波效果會(huì)變差[24],尤其數(shù)據(jù)含脈沖性質(zhì)的噪聲影響時(shí),系統(tǒng)的魯棒性欠佳。

針對(duì)非高斯白噪聲系統(tǒng)的狀態(tài)估計(jì),常用方法有狀態(tài)擴(kuò)展、量測(cè)差分[25]。但對(duì)于重尾分布的噪聲,很難基于一步相關(guān)模型去描述,同時(shí)也需要比較大的計(jì)算量,很難在實(shí)時(shí)應(yīng)用。多模型方法[26]也是針對(duì)非高斯白噪聲系統(tǒng)的手段,但主要缺點(diǎn)是多個(gè)模型并行運(yùn)行,計(jì)算量較大。

最大相關(guān)熵準(zhǔn)則(maximum correntropy criterion,MCC)作為非高斯信號(hào)處理的最佳魯棒準(zhǔn)則,近年來(lái)已成功應(yīng)用于自適應(yīng)濾波領(lǐng)域[27-28]。文獻(xiàn)[29]結(jié)合最大相關(guān)熵準(zhǔn)則和不動(dòng)點(diǎn)迭代算法,提出了最大相關(guān)熵卡爾曼濾波器(maximum correntropy Kalman filter,MCKF)。理論和實(shí)踐均表明,該濾波器對(duì)非高斯噪聲具有良好的魯棒性。

文中結(jié)合“當(dāng)前”概率模型,總結(jié)算法流程概況如下:

步驟2根據(jù)狀態(tài)預(yù)測(cè)和方差預(yù)測(cè)獲得一步預(yù)報(bào)值，即

P(k|k-1)=

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

計(jì)算P(k|k-1)的Cholesky分解獲得Bp(k|k-1);計(jì)算觀測(cè)噪聲陣R的Cholesky分解獲得Br(k)。

(1)

其中

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(12)

步驟5比對(duì)當(dāng)前狀態(tài)估計(jì)量和步驟4迭代狀態(tài)量,若滿(mǎn)足(13)式的迭代收斂條件則終止迭代,進(jìn)行步驟6,否則t+1→t,繼續(xù)步驟4迭代計(jì)算。

(13)

步驟6更新最終方差陣,k+1→k,回到步驟2 進(jìn)行下一點(diǎn)計(jì)算。

4 結(jié)果分析

4.1 參數(shù)設(shè)置

文中方法使用參數(shù)設(shè)置及原則如下:

(1) 目標(biāo)機(jī)動(dòng)常數(shù)的倒數(shù)α。根據(jù)表2中jounce頻率分布范圍,α取值范圍為0.5～10,結(jié)果差別不大,根據(jù)頻譜特征最集中段,分析時(shí)使用α=3。

(2) 觀測(cè)噪聲R。根據(jù)表1中噪聲統(tǒng)計(jì)結(jié)果可設(shè)對(duì)應(yīng)的觀測(cè)噪聲量。

(3) jounce機(jī)動(dòng)最大值自適應(yīng)參數(shù)。詳見(jiàn)第2.2節(jié)。

(4) 相關(guān)熵參數(shù)設(shè)置。包括核函數(shù)帶寬σ和迭代收斂條件ε。文獻(xiàn)[29]中指出,對(duì)于非脈沖噪聲,σ選取決定其性能是否能超過(guò)普通卡爾曼濾波器。以下分析時(shí),綜合考慮效果和性能,選擇σ=3,ε=10-6。

4.2 實(shí)際效果

衡量濾波結(jié)果好壞的指標(biāo)通常有均方根誤差和平滑度,其中平滑度可以限制過(guò)擬合,但對(duì)于船搖數(shù)據(jù)平滑度指標(biāo)設(shè)置比較困難,預(yù)測(cè)誤差也可以反映模型準(zhǔn)確度,限制過(guò)擬合,故文中結(jié)合濾波均方根誤差和三步預(yù)測(cè)誤差值的均方根誤差綜合考量。

圖6以橫搖數(shù)據(jù)為例給出了濾波效果,表3給出了濾波和預(yù)報(bào)誤差統(tǒng)計(jì)結(jié)果。

圖6 橫搖數(shù)據(jù)濾波及預(yù)報(bào)效果圖Fig.6 Filtering and forecasting effect of rolling data

表3 濾波及三步預(yù)報(bào)誤差均方根Table 3 Mean square error of filter and three-step prediction

可以得出以下結(jié)論:

(1) 該算法能很好地跟蹤船搖數(shù)據(jù),濾波誤差中不存在明顯趨勢(shì)量;濾波后數(shù)據(jù)平滑度較高,能有效濾除船搖數(shù)據(jù)噪聲。

(2) 濾波誤差統(tǒng)計(jì)結(jié)果與表1分析結(jié)果一致性較好。

(3) 表3的三步預(yù)報(bào)誤差統(tǒng)計(jì)結(jié)果表明,該算法能實(shí)現(xiàn)船搖數(shù)據(jù)的準(zhǔn)確預(yù)報(bào),均方根誤差小于5″。

(4) 表3的統(tǒng)計(jì)結(jié)果表明,該方法能準(zhǔn)確預(yù)報(bào)船搖角速度及角加速度,角速度預(yù)報(bào)均方根誤差約為0.006(°)/s,角加速度預(yù)報(bào)均方根誤差小于0.4(°)/s2。

(5) 圖6(d)表明,該方法魯棒性較好,尤其針對(duì)脈沖噪聲數(shù)據(jù),不會(huì)出現(xiàn)普通卡爾曼濾波器的精度變差的現(xiàn)象。

5 結(jié) 論

本文通過(guò)分析航天測(cè)量船的船搖數(shù)據(jù)特性,引入機(jī)動(dòng)目標(biāo)跟蹤理論進(jìn)行船搖濾波及預(yù)報(bào),針對(duì)現(xiàn)有模型階數(shù)不夠高的問(wèn)題導(dǎo)出了四階的CS-Jounce模型,并給出一種機(jī)動(dòng)參數(shù)自適應(yīng)模型。針對(duì)船搖數(shù)據(jù)環(huán)境非高斯白噪聲的問(wèn)題,引入了最大相關(guān)熵卡爾曼濾波器。實(shí)測(cè)數(shù)據(jù)分析證明,該方法能很好地跟蹤船搖數(shù)據(jù),濾波效果較理想,預(yù)報(bào)精度較高,且能很好地估計(jì)并預(yù)報(bào)船搖角速度及角加速度。