牛德智， 陳長(zhǎng)興， 班 斐， 陳 芳， 王 卓， 陳 強(qiáng)

(1.空軍工程大學(xué) 理學(xué)院,陜西 西安 710051；2.空軍工程大學(xué) 科研部,陜西 西安 710051;3.空軍工程大學(xué) 訓(xùn)練部,陜西 西安 710051)

動(dòng)態(tài)系統(tǒng)中的目標(biāo)跟蹤就是通過(guò)一組目標(biāo)的觀測(cè)量對(duì)其狀態(tài)量進(jìn)行準(zhǔn)確估計(jì)和預(yù)測(cè)的過(guò)程[1-2]，其核心部分是濾波算法?？柭鼮V波(KF)方法能有效跟蹤線性運(yùn)動(dòng)目標(biāo)，對(duì)于非線性問(wèn)題，擴(kuò)展卡爾曼濾波(EKF)和無(wú)跡卡爾曼濾波(UKF)分別具有一階和二階精度，但其性能會(huì)受非線性程度高低影響,且基于高斯假設(shè)條件[1]。近年來(lái)出現(xiàn)的粒子濾波(PF)方法用粒子集來(lái)近似后驗(yàn)概率密度，具有良好的非線性非高斯特性，廣泛應(yīng)用于目標(biāo)跟蹤、定位和導(dǎo)航、狀態(tài)監(jiān)視與故障診斷、圖像處理等方面[2]。

為解決傳統(tǒng)PF方法重采樣過(guò)程中的粒子退化問(wèn)題，出現(xiàn)了優(yōu)化組合重采樣、分層轉(zhuǎn)移重采樣等方法[3]，但是這些大都存在復(fù)雜度高、不易實(shí)現(xiàn)的缺點(diǎn)。另外，實(shí)際中的觀測(cè)噪聲多為非高斯噪聲，噪聲特性表現(xiàn)為概率密度函數(shù)具有“厚尾”特征，但大多數(shù)濾波算法在討論對(duì)非高斯噪聲的適應(yīng)性時(shí)，都是用特定概率密度分布(如gamma分布、瑞利分布[4]等)的噪聲來(lái)代替，并沒(méi)有用具有真實(shí)特征的非高斯噪聲來(lái)分析驗(yàn)證。

基于此，本文首先給出非高斯噪聲的產(chǎn)生方法，并從克服基本PF方法粒子退化現(xiàn)象的角度出發(fā)，提出一種易于實(shí)現(xiàn)且性能較優(yōu)的變步長(zhǎng)重采樣粒子濾波方法用于目標(biāo)跟蹤。

1 目標(biāo)跟蹤的基本原理

目標(biāo)跟蹤是利用傳感器所獲得的量測(cè)數(shù)據(jù)對(duì)目標(biāo)的狀態(tài)進(jìn)行持續(xù)和優(yōu)化的估計(jì)與預(yù)測(cè)過(guò)程，其中目標(biāo)狀態(tài)包括運(yùn)動(dòng)參數(shù)和屬性信息。實(shí)現(xiàn)原理[1]如圖1所示。

圖1 目標(biāo)跟蹤基本原理圖

對(duì)于一般的目標(biāo)跟蹤問(wèn)題，可有如下的目標(biāo)運(yùn)動(dòng)和觀測(cè)方程[5]：

式中,f(·)和 h(·)分別表示狀態(tài)轉(zhuǎn)移函數(shù)和觀測(cè)函數(shù)，xk為系統(tǒng)的狀態(tài)向量，zk為觀測(cè)向量，wk和 vk分別為過(guò)程噪聲和量測(cè)噪聲。濾波的過(guò)程就是從觀測(cè)序列zk中估計(jì)狀態(tài)向量 xk。

2 非高斯噪聲的產(chǎn)生方法

在已有的文獻(xiàn)資料中指出，可以由兩種以上的概率密度函數(shù)加權(quán)組合得到非高斯噪聲的概率密度函數(shù)[6]，如采用高斯概率密度函數(shù)和拉普拉斯概率密度函數(shù)組合就是一種方式，但是就其如何仿真形成具有實(shí)值序列的非高斯噪聲并沒(méi)有給出系統(tǒng)的方法流程；或者將具有兩種函數(shù)特性的噪聲序列值進(jìn)行加權(quán)組合，以替代加權(quán)組合概率密度函數(shù)作為實(shí)際操作，但這種仿真過(guò)程較為粗糙。為更好地體現(xiàn)非高斯噪聲的“厚尾”特性，本文給出兩種非高斯噪聲的產(chǎn)生方法。

將非高斯噪聲的概率密度函數(shù)表示為：

其中,PG(x)為高斯高斯概率密度函數(shù)，其表達(dá)式為PG(x)，通常情況下取標(biāo)準(zhǔn)正態(tài)分布,令μ=；PL(x)為拉普拉斯概率密度函數(shù)，其具有“厚尾”特性，表達(dá)式為；ε為加權(quán)系數(shù)，且 ε+(1-ε)=1保證了滿(mǎn)足 FNG(+∞)=1的條件限制。為了便于后文表述及在Matlab仿真中的應(yīng)用實(shí)現(xiàn)，聲明如下定義函數(shù)：

則非高斯噪聲的概率分布函數(shù)為：

2.1 概率分布函數(shù)產(chǎn)生法

首先求解F1(x)：

式(7)中的第4個(gè)等號(hào)表示可以借助Matlab中的函數(shù)命令求解反常積分的數(shù)值。

(1)當(dāng) x＞0 時(shí)：

(2)當(dāng) x≤0 時(shí) ：

一般地，ε描述了“厚尾”特性的概率密度函數(shù)所占權(quán)值大小，b表征了拉普拉斯概率分布的方差特性，如可取一組參數(shù)為：ε=0.3，b=2。這樣，可以得到一種生成非高斯噪聲序列的方法步驟如下：

(1)隨機(jī)生成 0～1之間的均勻分布隨機(jī)數(shù) U(0,1)，并令r(i)=U(0,1)，i=1,2,…,N為非高斯噪聲序列的序號(hào)值。

(3)令 X=[X(1)，X(2),… ,X(i),… ,X(N)],則 序 列 X 就是仿真生成的非高斯噪聲序列。

該方法為非高斯噪聲序列的仿真實(shí)現(xiàn)提供了一種工程借鑒意義。但是考慮到求解無(wú)解析形式原函數(shù)的積分會(huì)產(chǎn)生一定的時(shí)間復(fù)雜度，下面考慮一種依托加權(quán)組合概率密度函數(shù)的時(shí)效性更高的非高斯噪聲序列仿真實(shí)現(xiàn)手段。

2.2 生成數(shù)據(jù)庫(kù)產(chǎn)生法

在一個(gè)大值區(qū)間 x∈[-a，a]按一定的步進(jìn)步長(zhǎng) △(如△=0.01)取遍區(qū)間上的所有x,然后用數(shù)值定積分的方法，按照式(6)逐一獲得一系列點(diǎn)集處的混合噪聲(非高斯噪聲)的概率分布函數(shù)，并繪制曲線。這樣,大值區(qū)間上的數(shù)值點(diǎn)集與其對(duì)應(yīng)的概率分布函數(shù)值的一一映射對(duì)之間便構(gòu)成了一個(gè)非高斯噪聲數(shù)據(jù)庫(kù)D,并記D=[xi;Fi]，i=1,2,…,N,其中 xi∈[-a，a]，F(xiàn)i由 xi根據(jù)式(6)計(jì)算得來(lái)。

在此基礎(chǔ)上，便可以形成由數(shù)據(jù)庫(kù)生成非高斯噪聲序列的方法步驟：

(1)隨機(jī)生成0～1之間的均勻分布隨機(jī)數(shù) U(0,1)，并令r(i)=U(0,1)。

(2)求解k=arg min|r(i)-Fi|，進(jìn)而從數(shù)據(jù)庫(kù)D中取出xk。

(3)X=[xk(1)，xk(2),…,xk(N)],則序 列 X 就 是 仿真 生成的非高斯噪聲序列。

3 變步長(zhǎng)重采樣粒子濾波方法

3.1 基于置信區(qū)間的變步長(zhǎng)重采樣原理

針對(duì)粒子退化[7]問(wèn)題，變步長(zhǎng)重采樣方法由此產(chǎn)生，但是在現(xiàn)有的研究中，對(duì)變步長(zhǎng)重采樣的實(shí)現(xiàn)過(guò)程和方法描述還存在模糊界限，且主觀性較大，步長(zhǎng)較小時(shí)收斂速度仍然緩慢，步長(zhǎng)較大時(shí)仍會(huì)發(fā)生退化問(wèn)題[8]。

這樣，本文給出一種基于置信區(qū)間的變步長(zhǎng)重采樣方法的實(shí)現(xiàn)過(guò)程。由假設(shè)檢驗(yàn)的原理可知，若隨機(jī)變量Y服從均值為、方差為的某分布，則由大數(shù)定律可知，變量 Y的置信度為 1-α的置信區(qū)間為-σYuα/2+σYuα/2]，若置信水平 α=0.05，則 uα/2≈3，這就是工程上常用的統(tǒng)計(jì)數(shù)據(jù)是否合理的“3σ準(zhǔn)則”。借鑒這一思想，在PF的重采樣過(guò)程中，對(duì)于小權(quán)值粒子的處理過(guò)程如下：

若 xk(j)為小權(quán)值粒子，xk(i)為大權(quán)值粒子，則對(duì) xk(j)的重新賦值過(guò)程為:

其中，σq表示由重要性采樣函數(shù)所描述的概率分布的標(biāo)準(zhǔn)差，rand表示任意0～1之間的隨機(jī)數(shù)。

3.2 新的PF方法步驟

采用以上的重采樣方案后，得到目標(biāo)跟蹤的新方法算法步驟為：

(1)初始化：取k=0，以目標(biāo)的初始運(yùn)動(dòng)狀態(tài)作為先驗(yàn)p(x0)分布，建立初始狀態(tài)樣本{x0(i),1/N}Ni=1，N為粒子數(shù)。

(2)系統(tǒng)狀態(tài)轉(zhuǎn)移：根據(jù)k-1時(shí)刻的樣本粒子{xk(i)}Ni=1，通過(guò)狀態(tài)轉(zhuǎn)移方程得到k時(shí)刻的樣本粒子{xk(i)}Ni=1。

(3)生成觀察樣本：

(4)計(jì)算信息殘差：

(5)權(quán)值計(jì)算并歸一化:按照選取的重要性函數(shù)q(x)，計(jì)算各個(gè)粒子的權(quán)值，歸一化， 得支撐點(diǎn)集

(6)目標(biāo)狀態(tài)后驗(yàn)輸出：由支撐點(diǎn)集計(jì)算后驗(yàn)濾波的狀態(tài)輸出，即:

(7)重采樣：根據(jù)有效數(shù)據(jù)長(zhǎng)度原則從{xk(i),ωk(i)}Ni=1重新生成N個(gè)樣本，具體過(guò)程如下：

①產(chǎn)生在0～1上均勻分布的隨機(jī)數(shù)uj；

③計(jì)算uj落入權(quán)值累加函數(shù)的區(qū)間情況，即滿(mǎn)足Cm≥uj＞Cm-1；

④按照3.1節(jié)中描述的變步長(zhǎng)法重采樣:這樣，得到新的粒子集{xk(j),1/N}Ni=1。

(8)返回到第(2)步，跟蹤系統(tǒng)接收下一幀。重復(fù)上述迭代過(guò)程，獲得目標(biāo)的運(yùn)動(dòng)軌跡。

4 仿真實(shí)驗(yàn)

4.1 非高斯噪聲序列的產(chǎn)生

這里，首先仿真由加權(quán)組合形成的非高斯噪聲概率密度函數(shù)曲線，結(jié)果如圖2所示，從中可以看出其具有較長(zhǎng)的“厚尾”特性(參數(shù)選取為 ε=0.3，b=2)。

圖2 3種概率密度函數(shù)曲線圖

其次，圖3的仿真曲線為按照式(4)采用數(shù)值定積分法求解概率分布函數(shù)后繪制的曲線。這樣，在該曲線對(duì)應(yīng)的非高斯噪聲數(shù)據(jù)庫(kù)中，即可仿真給出具有一定長(zhǎng)度的非高斯噪聲序列。

圖4所示為由加權(quán)組合概率密度函數(shù)生成非高斯噪聲數(shù)據(jù)庫(kù)的仿真實(shí)現(xiàn)方法模擬生成的長(zhǎng)度為500的非高斯噪聲序列?？梢钥闯?，該噪聲序列具有良好的尖峰特性。

4.2 非線性運(yùn)動(dòng)目標(biāo)的跟蹤實(shí)驗(yàn)

圖3 數(shù)值定積分法繪制的3種噪聲概率分布函數(shù)

圖4 仿真模擬的長(zhǎng)度為500的非高斯噪聲序列

考慮使用雙坐標(biāo)雷達(dá)對(duì)平面上的一個(gè)運(yùn)動(dòng)目標(biāo)進(jìn)行觀測(cè)。目標(biāo)的起始點(diǎn)為(2 000,10 000),目標(biāo)在t=0～400 s沿y軸作恒速直線運(yùn)動(dòng)，運(yùn)動(dòng)速度為-15 m/s;在 t=400～600 s向x軸作90°的慢轉(zhuǎn)彎，加速度在 x方向和y方向的分量均為0.075 m/s2，完成慢轉(zhuǎn)彎后加速度將降為0;從 t=610 s開(kāi)始做 90°的快轉(zhuǎn)彎，加速度為 0.3 m/s2;在660 s結(jié)束轉(zhuǎn)彎，加速度降至0。雷達(dá)掃描周期 T=2 s，x和y獨(dú)立地進(jìn)行觀測(cè)，觀測(cè)噪聲的標(biāo)準(zhǔn)差均為100 m。對(duì)該問(wèn)題采用本文方法進(jìn)行跟蹤性能仿真,得到的仿真結(jié)果如圖5、圖 6所示。

圖5所示為新PF方法和標(biāo)準(zhǔn)PF方法(重采樣采用系統(tǒng)重采樣)的跟蹤性能，可知新PF方法對(duì)該目標(biāo)的跟蹤具有更好的效果和精度，由圖6中的跟蹤位置誤差便可說(shuō)明這一點(diǎn)。

圖5 兩種方法的跟蹤軌跡

圖6 兩種方法在x方向的跟蹤位置誤差

綜上可知，所提出的基于變步長(zhǎng)重采樣的新PF方法相比于標(biāo)準(zhǔn)PF方法，具有更好的跟蹤效果及精度，說(shuō)明其對(duì)于該類(lèi)雷達(dá)目標(biāo)跟蹤問(wèn)題具有一定的適應(yīng)性。

本文主要研究了非高斯噪聲條件下非線性運(yùn)動(dòng)目標(biāo)的跟蹤方法。首先給出了非高斯噪聲的仿真生成方法，進(jìn)而提出了一種基于置信區(qū)間的變步長(zhǎng)重采樣PF濾波方法。仿真實(shí)驗(yàn)結(jié)果表明，該方法可以有效應(yīng)用于非高斯噪聲條件下的非線性運(yùn)動(dòng)目標(biāo)的跟蹤，說(shuō)明了方法的有效性。

[1]Liu Changyun,Shui Penglang,Li Song.Unscented extended Kalman filter for target tracking[J].Journal of Systems Engineering and Electronics,2011,22(2):188-192.

[2]顧裕豐,劉國(guó)棟.一種復(fù)雜場(chǎng)景下的運(yùn)動(dòng)目標(biāo)跟蹤算法[J].電子技術(shù)應(yīng)用，2012,38(1):122-125.

[3]田雋,錢(qián)建生,李世銀.基于分層轉(zhuǎn)移的粒子濾波 MCMC重采樣算法[J].控制與決策,2011,26(8):1253-1258.

[4]孫景樂(lè),唐林波,趙保軍,等.基于瑞利分布的粒子濾波跟蹤算法[J].電子與信息學(xué)報(bào),2013,35(4):763-769.

[5]牛德智，陳長(zhǎng)興，王博,等.基于 FastICA的雷達(dá)弱目標(biāo)檢測(cè)前跟蹤方法研究[J].計(jì)算機(jī)應(yīng)用研究，2012,29(8):3120-3123.

[6]鄒國(guó)輝,敬忠良,胡洪濤.基于優(yōu)化組合重采樣的粒子濾波算法[J].上海交通大學(xué)學(xué)報(bào),2006,40(7):1135-1139.

[7]陳志敏,薄煜明,吳盤(pán)龍.一種新的混合智能粒子濾波算法在雷達(dá)機(jī)動(dòng)目標(biāo)跟蹤中的應(yīng)用[J].兵工學(xué)報(bào),2012,33(1):83-88.

[8]李秀智,劉紅云,居鶴華.粒子濾波趨優(yōu)重采樣算法及仿真研究[J].計(jì)算機(jī)仿真,2011,28(2):221-224.