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

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

0 引 言

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

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

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

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

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

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

1.1 時域分析

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

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

1.2 頻域分析

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

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

1.3 噪聲特性分析

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

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

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

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

表1 船搖數(shù)據(jù)噪聲均方差統(tǒng)計結(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ù)高階量信息可以通過差分的方法獲得。圖4橫搖為例,繪制了濾波后數(shù)據(jù)的一階至四階差分曲線圖。

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

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

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

表2 船搖數(shù)據(jù)jounce統(tǒng)計結(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ū)動過程,同時必須考慮個別時間段jounce振蕩加劇的現(xiàn)象。

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

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

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

2.1 CS-Jounce模型

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

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

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

量測方程為

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

當(dāng)采樣周期為T時,狀態(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足夠小時,矩陣退化為

其中

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足夠小時,Q(k)退化為

Q(k)=

jounce的預(yù)測方差為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(13)

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

4 結(jié)果分析

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

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

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

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

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

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

4.2 實(shí)際效果

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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