葉利民
(海軍工程大學研究生院,湖北武漢 430033)
隨著反艦導彈性能和突防技術的不斷提高,水面艦艇攔截反艦導彈變得更為困難。目前世界上許多國家均開始研發速度達2 Ma以上的超音速反艦導彈,而且攻擊模式也趨于多樣化,如水平機動、掠海飛行、躍升俯沖攻擊等,其中,末端超音速、躍升俯沖攻擊模式是現代導彈極具殺傷力的突防方案,而當前的目標跟蹤算法由于預置模型的局限性和目標運動的不確定性,對諸如上述高速強機動攻擊模式的跟蹤性能大大降低,無法滿足艦艇反導作戰的要求[1-2]。
目標自動跟蹤系統顯著提高了武器系統在戰場上的作戰能力與適應能力,它所采用的機動目標模型應是多模態的。機動目標參數辨識模型正是一種多模態模型,它是利用自適應濾波的思想,通過狀態濾波估計的輸出特性,構建信息反饋通道,動態地辨識機動目標的機動參數,但其對于非線性較強的系統辨識能力較差。粒子濾波算法是非常適合于非線性、非高斯系統的一種濾波算法。粒子濾波算法于20世紀50年代提出,60年代得到進一步發展,但一直受到粒子數匱乏和計算機處理能力的限制。直到1993年,Gordon等提出了一種新的基于Bootstrap的粒子濾波算法[3],才真正奠定了粒子濾波算法的基礎。粒子濾波的思想本質是用1組帶有權重值的隨機采樣點的離散和來計算難于處理的積分形式的后驗概率密度函數,并基于這些采樣點和權重值來計算狀態值。根據大數定理,當隨機采樣點足夠多時,粒子濾波算法可以基本不失真的表示狀態的后驗概率密度函數。粒子濾波算法也包括預測階段和更新階段,與卡爾曼濾波算法不同的是,它不需要依賴于系統方程和測量方程的線性程度和狀態的分布,克服了以往基于線性高斯濾波的缺點,其主要優點包括:
1)可以表示任意的概率密度函數;
2)對狀態空間的可能區域進行自適應的關注;
3)可以處理非高斯的系統噪聲;
4)濾波框架可以實現對機動目標的跟蹤。
本文通過將以往參數辨識過程中的卡爾曼濾波器換為粒子濾波器,發揮粒子濾波可以對非線性非高斯噪聲系統進行濾波的特點,結合參數辨識的遞推最小二乘算法,提出了一種機動目標自適應跟蹤算法。仿真實驗說明了該算法的有效性。
在連續的時間域中,以三階線性微分方程結構確定待估計的機動目標運動模型的結構[4]為:

其中:a和b為待估計的模型參數;W(t)是均值為0的白噪聲。
由于要通過模型參數a和b的估計來實現模型的估計,它的狀態方程為:

由參數估計模型的表達式(2)可以看出,參數a和b綜合決定了目標加速度的變化率。b是變化率中的常數,單位是m/s3,a是相對于加速度的機動系數,一般情況下是1個小于0的負數,單位是1/s。
依據參數a和b的不同組合,將使參數估計模型的結構可以覆蓋機動目標可能的多種運動狀態。
上述參數估計模型對目標的多種運動模態的廣泛適應能力,使得它不僅具有對采樣周期變化的適應能力,而且在它的有效動態遞推估計算法中通過參數a和b的自動校正,還能自動適應機動目標運動模態的變化,不需要像IMM算法那樣,需要對機動目標多種運動模型之間的轉換規律進行人為的假定和估計。上述模型的特性為自適應濾波的實現提供了基本條件,利用目標運動信息估計得到參數a和b,獲得下次濾波的模型,原理如圖1所示。
依據參數估計模型(2),以求解微分方程的歐拉算法的二步格式和中心差商的數值微分算法為基礎,實現微分方程的動態估計。在一維的情況下,將待估計的參數估計模型的結構表示為:

圖1 自適應濾波原理框圖Fig.1The principle of adaptive filter

其中,e(t)為模型的殘差。
模型的狀態方程為:

根據求解微分方程的二步格式的歐拉算法,選定步長為h,對狀態方程式(4)進行線性變換,

其中:k為步長為h的時間序號。
利用差商數值微分公式有

將式(6)代入差分方程組(5),由下向上迭代化簡得:

令x(k)=x1(k),得到以待估計參數a和b為變量的線性方程:

線性方程式(8)中各已知量之間的變換關系式為:

當系統動態觀測到目標位移量序列值x(-3),x(-2),x(-1),x(0),…,x(m+1),m≥2,則可按式(8)順序構造出如下線性方程組:

極小化指標函數J,即可得到目標“參數估計模型”參數θ的最小二乘解

只要UTU非奇異,就有惟一解,這在通常的動態觀測中是可以滿足的。因此,利用線性變換式(9),就可實現“參數估計模型”的動態估計。而且,模型參數θ的估計,還可采用遞推算法,這正是自適應濾波所必需的。
對于用m個方程所得的解,記為


基于貝葉斯濾波思想,粒子濾波(Particle Filters)通過帶有歸一化權值的粒子集來近似表示后驗概率密度。當粒子采樣數量足夠大時,能準確地表達后驗密度分布,此時的粒子濾波算法接近貝葉斯最優估計。每個粒子的位置和權重反映了狀態空間在該區間的密度[5-7]。設采樣N個粒子粒子濾波算法原理描述如下式:


在通常情況下,建議將密度函數分解為:

根據貝葉斯公式,后驗概率密度可表示為:

需要說明的是,與上面討論的遞推貝葉斯估計的假設相同,這里的建議密度函數同樣服從一階馬爾科夫過程,且已知觀測獨立于指定狀態。根據貝葉斯估計原理,粒子濾波的過程可分為以下2個階段:
1)預測階段
預測階段是先驗概率在狀態空間進行搜索的過程。從過程噪聲ωk-1中采樣N個點,使用這些點通過狀態方程形成新的粒子分布,其計算公式為:

2)更新階段
當得到k時刻的觀測值,將式(20)和式(21)代入權重計算式(19),可以得到粒子權重

通常需要對上面計算的權重進行歸一化得到權重


上述方程式構成了遞推的PF算法的基本過程。通過上述推導過程可以看出,粒子濾波不受線性高斯條件的約束,易于非線性非高斯跟蹤算法的實現。
通過將參數估計的遞推最小二乘濾波算法和粒子濾波算法結合起來,可以實現自適應估計濾波,由濾波輸出結果估計得到參數a和b,由參數a和b確定粒子濾波算法中目標模型的狀態轉移矩陣,如此循環計算,算法流程如圖2所示。
利用Matlab軟件,采用蒙特卡洛法進行仿真實驗。

圖2 基于參數估計模型的濾波算法實現流程圖Fig.2The flow chart of realizing the filtering algorithm based on parameter identification model
設目標在做角速度恒定的轉彎運動:初始斜距600 m;距我艦初始方位角0.785 rad;變化率1 rad/s;速度為850 m/s;雷達采樣率為50 Hz;采樣持續時間10 s;觀測距離隨機誤差5 m;觀測方位角隨機誤差0.9 mrad。
分別用CT跟蹤模型以及基于粒子濾波算法和遞推最小二乘算法的目標參數估計模型對目標進行跟蹤,仿真結果如圖3~圖6所示。
由仿真結果可以發現,本文提出的算法對機動目標的跟蹤效果要好于用單一模型進行跟蹤的效果。本文所提模型既可動態地進行識別,又具有對工作周期變化的適應能力,可以覆蓋機動目標可能的多種運動狀態,使自適應濾波問題可在正常的技術思路下進行。仿真結果表明,此模型在原理上是正確的,在計算上是可行的,濾波器能較快收斂,計算精度高,反應時間短,在解決機動目標跟蹤問題上,效果較好。

圖3 x軸濾波效果Fig.3The filter performance in x axis

本文針對反艦導彈末端機動問題,通過對參數自動辨識模型進行改進,引入對非線性非高斯噪聲系統有良好跟蹤效果的粒子濾波算法,提出了一種新的參數辨識模型。并在目標作勻速圓周運動的情況下進行了算法的跟蹤性能仿真試驗,試驗結果證實了算法的有效性,具有一定工程實際意義。
[1]LI X R,JILKOV V P.A survey of maneuvering target tracking partⅠ:dynamic models[J].IEEE Trans.Aerospace and Electronic Systems(S0018-9251),2003,39(4):1333-1363.
[2]LI X R,JILKOV V P.A Survey of maneuvering target tracking partⅡ:ballistic target Models[C]//Proc.2001 SPIE Conf.on Signal and Data Processing of Small Targets,San Diego,CA,USA.USA:SPIE Press,2001,4473:559-581.
[3]GORDON N J,SALMOND D J,SMITH A F M.Novel approachtononlinear/non-GaussianBayesianstate estimation[J].Proceeding of the IEEE,1993,140(2):107-113.
[4]王志賢.最優狀態估計與系統辨識[M].西安:西北工業大學出版社,2004.181-219.
[5]CRISAN D,DOUCET A.A survey of convergence results on particle filtering methods for practitioners[J].IEEE Transaction on Signal Processing,2002,50(2):736-746.
[6]KALMAN H S.Filtering and neural networks[M].New York:John Wiley and Sons,2001.
[7]LIU J S,Chen R.Sequential monte carlo methods for dynamical systems[J].Journal of the American Statistical Association,1998,93(5):1032-1044.