郭夢(mèng)春，駱吉安

(杭州電子科技大學(xué)自動(dòng)化學(xué)院，浙江 杭州 310018)

0 引 言

基于時(shí)延和多普勒頻移的被動(dòng)輻射源狀態(tài)估計(jì)廣泛應(yīng)用于雷達(dá)、聲納、衛(wèi)星定位和導(dǎo)航等領(lǐng)域，接收站利用自身的位置以及接收到的輻射源信號(hào)實(shí)時(shí)解算出目標(biāo)的位置和速度等狀態(tài)參數(shù)[1]。傳統(tǒng)兩步跟蹤算法首先估計(jì)中間參數(shù)，然后解算中間參數(shù)構(gòu)成的觀測(cè)方程從而獲得目標(biāo)的跟蹤軌跡，如到達(dá)角(Angle of Arrival, AOA)[2-3]、到達(dá)時(shí)間(Time of Arrival, TOA)[4-5]、到達(dá)時(shí)間差(Time Difference of Arrival, TDOA)[6-7]、多普勒頻移[8-9]和到達(dá)頻率差(Frequency Difference of Arrival, FDOA)[10]等是常見(jiàn)的中間參數(shù)估計(jì)方法在兩步法中的應(yīng)用。但是，兩步法在參數(shù)估計(jì)階段忽略了觀測(cè)站之間的聯(lián)系性，損失了部分位置信息，同時(shí)引入的誤差不斷疊加累積，因此，兩步跟蹤算法是次優(yōu)的。

目標(biāo)直接跟蹤算法省去中間參數(shù)估計(jì)步驟，直接使用固定或移動(dòng)接收器截獲的信號(hào)波形來(lái)跟蹤移動(dòng)發(fā)射器。文獻(xiàn)[11]采用多普勒頻率測(cè)量對(duì)勻速運(yùn)動(dòng)目標(biāo)進(jìn)行直接定位，但該方法僅適用于基于多普勒頻移的定位機(jī)制，無(wú)法推廣至AOA，TDOA等其他直接定位方法。文獻(xiàn)[12]提出一種基于時(shí)延和多普勒頻移的直接跟蹤粒子濾波算法，在處理時(shí)延信息時(shí)采用向下取整的方式代替真實(shí)的時(shí)延信息，舍棄了小數(shù)部分并引入了量化誤差，同時(shí)粒子濾波器的重要性重采樣降低了粒子群的多樣性，引入了粒子的貧化問(wèn)題，影響其跟蹤精度。為此，文獻(xiàn)[13]在處理時(shí)延信息時(shí)，采用離散傅里葉變換和反變換，沒(méi)有舍棄其小數(shù)部分，避免了量化誤差的引入，提高了跟蹤精度，但忽略了信號(hào)傳播過(guò)程中衰減系數(shù)的估計(jì)，使用的粒子濾波器仍存在粒子貧化問(wèn)題。為了避免粒子多樣性減少問(wèn)題，本文采用馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo, MCMC)粒子直接跟蹤算法來(lái)實(shí)現(xiàn)運(yùn)動(dòng)目標(biāo)的位置、速度以及高度的聯(lián)合估計(jì)，在粒子濾波迭代中引入典型的MCMC算法—M-H抽樣算法，有效應(yīng)對(duì)了粒子貧化問(wèn)題，提高了跟蹤精度。

1 問(wèn)題描述

rk+1=Frk+wk

(1)

式中，F(xiàn)為目標(biāo)的狀態(tài)轉(zhuǎn)移矩陣，wk為均值為零協(xié)方差矩陣為Qw的高斯白噪聲過(guò)程，F(xiàn)和Qw分別為：

F=diag{[F′,F′,F′]}

(2)

(3)

式中，q為系統(tǒng)定義的參數(shù)，F(xiàn)′和Q′為中間變量，分別為：

(4)

(5)

式中，Tr為采樣時(shí)間間隔，假設(shè)狀態(tài)模型的條件概率密度函數(shù)p(rk+1|rk)服從高斯分布，有：

p(rk+1|rk)～N(Frk,Qw)

(6)

最后假設(shè)初始狀態(tài)r1具有已知的先驗(yàn)高斯分布為：

r1～N(r0,Qr0)

(7)

式中，r0為給定的初始狀態(tài)向量，Qr0為給定的對(duì)角矩陣。

假設(shè)在空間內(nèi)有M個(gè)固定的接收站，它們?cè)跁r(shí)間和頻率上是同步的。接收站在K個(gè)觀測(cè)間隔內(nèi)對(duì)輻射源的信號(hào)進(jìn)行截獲，um,k=[xm,k,ym,k,zm,k]T表示第m個(gè)接收站在第k觀測(cè)時(shí)刻的位置坐標(biāo)，設(shè)采樣點(diǎn)數(shù)為nk，則第m個(gè)接收站在第k時(shí)刻接收的采樣信號(hào)為：

(8)

式中，j表示虛部，Ts為采樣周期，bm,k為未知的路徑衰減系數(shù)，sk(t)為第k個(gè)觀測(cè)間隔目標(biāo)的復(fù)數(shù)信號(hào)，nm,k(t)是零均值的廣義平穩(wěn)的復(fù)高斯白噪聲，其協(xié)方差矩陣為Γk=σn2I，τm,k為信號(hào)到達(dá)第m個(gè)接收站的時(shí)延，fm,k為第m個(gè)接收站與目標(biāo)相對(duì)位移引起的多普勒頻移，分別為：

(9)

(10)

zm,k=bm,kAfm,kGHAτm,kGsk+nm,k=bm,kHm,ksk+nm,k

(11)

式中，G為離散DFT矩陣，GH為矩陣G的共軛轉(zhuǎn)置，Afm,k為多普勒頻移矩陣，Aτm,k為時(shí)延矩陣，其表達(dá)式分別為：

(12)

Afm,k=diag{exp(j2πTsfm,knk)}

(13)

(14)

聯(lián)立所有觀測(cè)方程有：

(15)

在發(fā)射信號(hào)和路徑傳播衰減系數(shù)未知的情況下，接收站直接截獲運(yùn)動(dòng)目標(biāo)輻射的無(wú)線電信號(hào)，建立關(guān)于目標(biāo)狀態(tài)的觀測(cè)方程。本文研究的目標(biāo)是通過(guò)已知的非線性觀測(cè)信息zk對(duì)運(yùn)動(dòng)目標(biāo)x軸方向位置rx,k，y軸方向位置ry,k，z軸方向高度rz,k以及目標(biāo)的速度進(jìn)行聯(lián)合精確估計(jì)。標(biāo)準(zhǔn)粒子濾波器在重采樣階段會(huì)導(dǎo)致粒子群出現(xiàn)單一化，影響跟蹤性能，為了提高跟蹤精度，本文在粒子濾波迭代中引入MCMC算法。

2 基于MCMC粒子直接跟蹤算法

直接跟蹤粒子濾波算法有一個(gè)局限性是粒子退化問(wèn)題。重采樣是抑制粒子退化的有效方法，通過(guò)在重要性采樣階段實(shí)施重采樣，淘汰權(quán)重低的粒子，集中復(fù)制權(quán)重高的粒子，從而達(dá)到抑制粒子退化目的。但是少數(shù)權(quán)重較高的粒子被大量復(fù)制從而降低了粒子群的多樣性[15]，導(dǎo)致采樣樣本枯竭，影響跟蹤性能。

為了避免粒子貧化問(wèn)題，本文在粒子濾波器順序迭代中引用MCMC算法估計(jì)目標(biāo)的后驗(yàn)概率分布。MCMC算法的核心思想是通過(guò)建立一個(gè)平穩(wěn)分布為p(rk|z1:k)的馬氏鏈得到近似服從p(rk|z1:k)的隨機(jī)樣本，然后采用這些樣本做各種統(tǒng)計(jì)，MCMC算法構(gòu)成的馬氏鏈具有很好的收斂性。由于直接跟蹤的觀測(cè)模型是高度非線性的，其似然函數(shù)可以近似表示為：

(16)

(17)

式中，Φk為M×M的矩陣，其表達(dá)式為：

Φk=(Uk)HUk

(18)

Uk=[(H1,k)Hz1,k,…,(HM,k)HzM,k]

(19)

MCMC粒子算法最重要部分是構(gòu)造轉(zhuǎn)移核，常見(jiàn)的構(gòu)造轉(zhuǎn)移核抽樣方法有M-H抽樣[16]和吉布斯抽樣(Gibbs Sampling)，基于M-H抽樣的MCMC粒子直接跟蹤算法主要步驟如下。

(2)從區(qū)間(0,1)中按照均勻概率分布抽樣得隨機(jī)值u，即u～U(0,1)。

(4)計(jì)算拒絕-接收概率

(20)

上述步驟中，步驟3中提議分布的選擇是最重要部分，理想的提議分布應(yīng)獨(dú)立于以p(rk|z1:k)為平穩(wěn)分布的馬氏鏈，其分布組成如下：

(21)

式中，Nb為MCMC算法預(yù)燒期的迭代次數(shù)，N為算法迭代保持?jǐn)?shù)(平穩(wěn)鏈的長(zhǎng)度)，Nb+N為算法迭代總次數(shù)。式(21)提議分布的采樣過(guò)程分為以下2個(gè)步驟。

(22)

(23)

上述采樣過(guò)程是理想的采樣方法，但是在大多數(shù)場(chǎng)合從p(rk|zk,rk-1)采樣是很困難的。常見(jiàn)的是選擇基于先驗(yàn)概率分布作為目標(biāo)的提議分布，最簡(jiǎn)單的替代方法是選擇先驗(yàn)分布和上一時(shí)刻的近似后驗(yàn)分布概率的組合作為目標(biāo)的提議分布，該提議分布表示為：

(24)

式(24)的采樣過(guò)程分為如下2個(gè)步驟進(jìn)行：

(25)

(26)

使用式(24)作為目標(biāo)的提議分布，基于M-H抽樣的MCMC粒子直接跟蹤算法的拒接-接收概率僅由似然函數(shù)比給出，式(20)改寫(xiě)為：

(27)

步驟4中的拒絕-接收概率ρ改成由式(27)定義。MCMC粒子直接跟蹤算法提議分布的選擇是重要步驟，理想的提議分布是從后驗(yàn)概率分布中直接采樣，但這種采樣方法非常困難。本文將上一時(shí)刻的后驗(yàn)分布與當(dāng)前時(shí)刻的先驗(yàn)分布組合起來(lái)作為目標(biāo)的提議分布，先從上一時(shí)刻的后驗(yàn)分布中采樣，再?gòu)漠?dāng)前時(shí)刻的先驗(yàn)分布中生成采樣點(diǎn)，有效解決了MCMC算法的采樣困難問(wèn)題。

3 仿真結(jié)果與分析

通過(guò)在三維空間內(nèi)部署多個(gè)接收站對(duì)移動(dòng)發(fā)射器的直接跟蹤來(lái)估計(jì)所提算法的性能。假設(shè)接收站個(gè)數(shù)M=4，接收站靜止且近似線性均勻排列，位置坐標(biāo)分別為：

(28)

(29)

圖1 位置跟蹤過(guò)程

設(shè)置馬爾科夫鏈總長(zhǎng)為500，得到三維跟蹤效果如圖1所示，運(yùn)動(dòng)方向從右向左。從圖1可以看出，MCMC粒子直接跟蹤算法經(jīng)過(guò)短暫的調(diào)整之后以一定的誤差穩(wěn)定地跟蹤目標(biāo)，因?yàn)镸CMC粒子直接跟蹤算法構(gòu)造了一個(gè)服從目標(biāo)后驗(yàn)概率密度函數(shù)p(rk|z1:k)的馬氏鏈，通過(guò)對(duì)馬氏鏈樣本作平均來(lái)實(shí)時(shí)估計(jì)目標(biāo)狀態(tài)。

為了進(jìn)一步說(shuō)明MCMC粒子直接跟蹤算法的跟蹤性能，將其與標(biāo)準(zhǔn)粒子濾波直接跟蹤算法[11]進(jìn)行比較。目標(biāo)位置的均方根誤差(Root Mean Square Error, RMSE)和速度RMSE分別為：

(30)

圖2 不同算法跟蹤誤差曲線圖

4 結(jié)束語(yǔ)

針對(duì)運(yùn)動(dòng)目標(biāo)直接跟蹤粒子濾波器存在粒子貧化問(wèn)題，本文提出一種基于時(shí)延和多普勒頻移的MCMC粒子直接跟蹤算法。在發(fā)射信號(hào)及路徑衰減系數(shù)均未知的情況下，使用時(shí)間延遲和多普勒頻移原理建立直接跟蹤觀測(cè)模型，并在粒子濾波器順序迭代中引入MCMC經(jīng)典算法—M-H抽樣算法，選擇滿足拒絕-接收率的采樣樣本作為新的粒子群，避免了傳統(tǒng)粒子濾波重采樣階段采樣枯竭問(wèn)題，進(jìn)一步提高了目標(biāo)跟蹤性能。但是，本文算法主要是在提高精度方面的突破，下一步將對(duì)復(fù)雜度較小的算法展開(kāi)研究。