基于動態(tài)規(guī)劃和雙變量容積卡爾曼濾波的空基外輻射源定位算法

盧雨，衣曉

(海軍航空大學(xué)，山東 煙臺 264001)

0 引言

空基外輻射源定位(APCL)是指利用空基外輻射源對目標(biāo)進(jìn)行定位跟蹤的一種技術(shù)[1-6]。區(qū)別于地基外輻射源如數(shù)字廣播電視[2-3]、移動手機(jī)信號[4]等，以及天基外輻射源如通信衛(wèi)星[5]、全球定位導(dǎo)航信號[6]等，空基外輻射源信號則是指由空中移動平臺所輻射的電磁信號。這種外輻射源信號不具備地基、天基兩種信號所擁有的穩(wěn)定性和可靠性，而且外輻射源信號是非合作的，空基外輻射源平臺的運(yùn)動狀態(tài)未知且難測，這就極大地影響了系統(tǒng)對目標(biāo)的定位跟蹤精度。因此，如何在復(fù)雜的雜波環(huán)境下提取目標(biāo)狀態(tài)，并降低外輻射源不確定性對系統(tǒng)定位跟蹤精度的影響是提升APCL系統(tǒng)實用效能的重要環(huán)節(jié)。

針對無源定位系統(tǒng)中存在傳感器狀態(tài)不確定性的影響，文獻(xiàn)[7]提出一種迭代最大似然概率數(shù)據(jù)關(guān)聯(lián)(IML-PDA)算法，該算法采用迭代預(yù)測-更新框架，有效提高了目標(biāo)的定位精度。針對外輻射源定位系統(tǒng)中存在外輻射源狀態(tài)不確定的影響，文獻(xiàn)[8]給出了全面詳細(xì)的應(yīng)對方案，該方案采用聯(lián)合狀態(tài)估計的方法，提出一種遞推的航跡起始/保持算法，將外輻射源狀態(tài)和目標(biāo)狀態(tài)同時進(jìn)行更新，當(dāng)目標(biāo)和外輻射狀態(tài)相互獨(dú)立時，可以取得較好的效果。實際應(yīng)用背景下，過程和量測噪聲并非服從高斯分布且統(tǒng)計特性不精確，針對這種情況，文獻(xiàn)[9]提出一種基于最大熵方法的變分貝葉斯自適應(yīng)卡爾曼濾波，進(jìn)而提高了系統(tǒng)的濾波精度和穩(wěn)定性。當(dāng)環(huán)境中存在強(qiáng)雜波干擾時，則需要考慮目標(biāo)是否存在，如何提取出目標(biāo)狀態(tài)，文獻(xiàn)[10]提出了一種基于序貫蒙特卡洛與概率假設(shè)密度濾波相結(jié)合的目標(biāo)跟蹤方法，在主動分布式聲納系統(tǒng)的目標(biāo)定位應(yīng)用中取得了較好的效果，但是該方法在被動跟蹤場景中的應(yīng)用效果還有待研究。文獻(xiàn)[11]基于極值理論確定航跡回溯的閾值，利用誤差傳播理論構(gòu)建動態(tài)規(guī)劃(DP)的狀態(tài)轉(zhuǎn)移范圍，并通過幅值積累的方式提取出了目標(biāo)的狀態(tài)，但是采用的狀態(tài)擴(kuò)維(SA)方法計算量較大，且沒有考慮目標(biāo)和外輻射源狀態(tài)的相互獨(dú)立性對目標(biāo)和外輻射源狀態(tài)估計造成的影響。針對APCL系統(tǒng)在實際應(yīng)用中可能出現(xiàn)的不可觀測問題，文獻(xiàn)[12]基于條件數(shù)定義的可觀測度對其進(jìn)行了詳細(xì)的分析。

考慮到APCL系統(tǒng)的應(yīng)用環(huán)境相對復(fù)雜，目標(biāo)回波功率相對微弱，相比于外輻射源的直達(dá)波而言，容易出現(xiàn)被雜波淹沒的問題。本文研究的主要問題就是如何提取雜波中的目標(biāo)量測，并降低外輻射源狀態(tài)不確定性的影響，實現(xiàn)對目標(biāo)和外輻射源的狀態(tài)估計。針對雜波環(huán)境下的空基外輻射源定位問題，本文提出一種基于DP和雙變量容積卡爾曼濾波(DP-BVCKF)的定位跟蹤算法。首先基于誤差傳播理論，確定出目標(biāo)的狀態(tài)轉(zhuǎn)移范圍；其次，選取回波信號幅度值作為 DP的值函數(shù)，在每個候選測量值的狀態(tài)轉(zhuǎn)移范圍中尋找相鄰時刻之間的最大幅度值進(jìn)行遞推積累；再次，基于雜波的統(tǒng)計特性和極值理論來解算DP的目標(biāo)提取閾值，逆向提取來自目標(biāo)的測量值序列[11]；然后利用雙變量容積卡爾曼濾波(BVCKF)算法同時對目標(biāo)和外輻射源進(jìn)行狀態(tài)估計，從原理上降低外輻射源狀態(tài)不確定性對目標(biāo)狀態(tài)估計造成的影響。

1 問題描述

Xk+1=FkXk+Gkuk+Vk，

(1)

強(qiáng)雜波環(huán)境下的APCL系統(tǒng)測量值將會包含源于目標(biāo)、外輻射源以及雜波的信息。假設(shè)k時刻觀測站接收到Nk個量測數(shù)據(jù)，其數(shù)據(jù)集合為

{zk(i)=

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

幅度值A(chǔ)k(i)的測量方程為

(10)

(11)

2 DP-BVCKF算法描述

當(dāng)目標(biāo)回波較弱時，難以區(qū)分雜波測量值和目標(biāo)測量值，再加上外輻射源狀態(tài)不確定性的影響，直接采用一般的濾波算法很容易形成虛假目標(biāo)航跡。針對APCL系統(tǒng)的上述問題，本文提出DP-BVCKF算法，該算法主要分兩個步驟：一是對一組已獲取的含雜波量測值，進(jìn)行DP[16]處理，逆向提取目標(biāo)的量測序列；二是采取BVCKF算法，同時對目標(biāo)和外輻射源進(jìn)行狀態(tài)估計，降低外輻射源狀態(tài)不確定性對目標(biāo)狀態(tài)估計的影響。具體算法流程如圖2所示。圖2中k+1|k、Px,k+1|k分別為目標(biāo)的狀態(tài)向量和協(xié)方差，k+1|k、Py,k+1|k分別為外輻射源的狀態(tài)向量和協(xié)方差，k+1|k為量測的一步預(yù)測值，h為量測轉(zhuǎn)移矩陣，Kx、Ky分別為目標(biāo)和外輻射源的卡爾曼濾波增益。

圖2 DP-BVCKF算法流程圖

2.1 構(gòu)建DP的狀態(tài)轉(zhuǎn)移范圍

在APCL系統(tǒng)的非線性測量空間中，當(dāng)前時刻的目標(biāo)在下一時刻可能出現(xiàn)的狀態(tài)集合記為目標(biāo)的狀態(tài)轉(zhuǎn)移范圍記為Ω(zk+1)?？紤]到外輻射源不確定性的影響，Ω(zk+1)的大小主要取決于θk(i)、φk(i)以及rk(i)，其范圍參數(shù)可以由目標(biāo)方位角和到達(dá)時間差的變化量來表示，即Ω(zk+1)=[Δθk,Δrk]。因此，根據(jù)APCL系統(tǒng)的非線性測量信息可以構(gòu)建目標(biāo)的狀態(tài)轉(zhuǎn)移范圍[11]為

(Δθk)2=(Δxk)2(?θk/?xk)2|xk=α(zk(i))+

(Δyk)2(?θk/?yk)2|yk=β(zk(i))，

(12)

(Δrk)2=(Δxk)2(?rk/?xk)2|xk=α(zk(i))+

(Δyk)2(?rk/?yk)2|yk=β(zk(i))+

(Δxt,k)2(?rk/?xt,k)2|xt,k=t,k+

(Δyt,k)2(?rk/?yt,k)2|yt,k=t,k，

(13)

α(zk(i))=

(14)

β(zk(i))=

(15)

2.2 幅值積累和航跡逆向提取

從雜波中提取目標(biāo)量測數(shù)據(jù)的核心在于雜波的幅值大小取決于其統(tǒng)計特性，而目標(biāo)的幅值受多種因素影響，在利用DP方法選取幅值A(chǔ)k(i)作為值函數(shù)[16]對相鄰時刻中的量測數(shù)據(jù)進(jìn)行幅值積累時，目標(biāo)的狀態(tài)序列所積累的值函數(shù)往往是最大的。根據(jù)極值理論，可以求取一定檢測概率下的雜波幅度極值，并將其作為幅值積累的門限閾值，逆向提取目標(biāo)的狀態(tài)序列。記DP的值函數(shù)為I(zk(j))，j表示用于值函數(shù)累積的第j個測量值。目標(biāo)狀態(tài)序列逆向提取的具體步驟如下：

1)確定值函數(shù)積累的門限閾值。根據(jù)極值理論，雜波幅值的極大值服從Gumbel分布，其函數(shù)表示式為

(16)

式中：FA(·)為Gumbel分布函數(shù)；PFT為虛警率；a和b為FA(·)參數(shù)，可通過(17)式、(18)式[17]迭代估計得到：

(17)

(18)

2)幅值積累初始化。當(dāng)k=1時，遍歷當(dāng)前時刻的所有量測數(shù)據(jù)，令值函數(shù)I(z1(i))=A1(i)，記η0(z1(i))=0為目標(biāo)狀態(tài)序列。

3)幅值積累，閾值判斷，終止積累。當(dāng)k>1時，遍歷k時刻的所有量測數(shù)據(jù)，進(jìn)行幅值積累，過程為

(19)

(20)

式中：Ωk-1(zk(i))為第k-1時刻量測數(shù)據(jù)可能轉(zhuǎn)移到zk(i)的狀態(tài)轉(zhuǎn)移范圍；ηk-1(zk(i))為轉(zhuǎn)移到zk(i)的第k-1時刻測量值。再根據(jù)步驟1求得的門限閾值S，提取出滿足判定條件的所有狀態(tài)序列集合M(i)，判定條件為

{M(i)}={z1∶M(i)∶I(zM(i))>S},

(21)

式中：z1∶M(i)表示滿足判定條件的M個狀態(tài)序列。

2.3 基于BVCKF算法的狀態(tài)估計

BVCKF算法的核心思想是將量測一步預(yù)測過程中對單變量的多維積分?jǐn)U展為對雙變量的二重多維積分，并在互協(xié)方差計算過程中對單變量的概率密度函數(shù)求解轉(zhuǎn)換為對雙變量聯(lián)合概率密度函數(shù)的邊緣概率密度函數(shù)的求解。由于非線性函數(shù)的多維積分不易求解[18]，該算法采用了球面-徑向容積準(zhǔn)則[19-20]來對非線性函數(shù)的多維積分進(jìn)行近似計算，既降低了計算量，又保證了較高的近似精度，具體算法流程描述如下：

(22)

(23)

(24)

式中：n為狀態(tài)向量的維數(shù)；ξi為初始容積點(diǎn)集；[1]i表示第i個初始容積點(diǎn)，初始容積點(diǎn)集[1]的具體形式為

[1]=

2)狀態(tài)一步預(yù)測。首先計算經(jīng)過系統(tǒng)狀態(tài)方程f(·)傳遞的容積點(diǎn)：

(25)

(26)

k+1|k(k+1|k)T+Γk，

(27)

式中：ωi為權(quán)重，ωi=1/(2n)。同理可得Y在k+1時刻的狀態(tài)一步預(yù)測k+1|k和一步預(yù)測協(xié)方差矩陣Py,k+1|k.

3)重新構(gòu)造容積點(diǎn)。根據(jù)球面-徑向容積準(zhǔn)則，當(dāng)?shù)玫綘顟B(tài)一步預(yù)測k+1|k及其對應(yīng)的協(xié)方差矩陣后，需要重新構(gòu)造容積點(diǎn)，以便更好地近似變量的概率密度函數(shù)：

(28)

(29)

4)雙變量聯(lián)合的量測一步預(yù)測。將對單一變量的多重積分?jǐn)U展為對雙變量的二重積分，再對二重積分進(jìn)行數(shù)值近似，得到雙變量聯(lián)合的量測一步預(yù)測及其協(xié)方差為

(30)

(31)

(32)

5)增益的計算。首先計算狀態(tài)變量X的互協(xié)方差：

(33)

(34)

再計算狀態(tài)變量的卡爾曼增益：

(35)

同理可得狀態(tài)變量Y對應(yīng)的卡爾曼增益Ky.

6)狀態(tài)更新。利用計算得到的卡爾曼增益和最新時刻的量測值就可以對狀態(tài)變量值和協(xié)方差矩陣進(jìn)行更新：

k+1|k+1=k+1|k+Kx(Zk+1-k+1)，

(36)

(37)

式中：Zk+1對應(yīng)3.2節(jié)中提取的M(i)中的目標(biāo)量測數(shù)據(jù)。同理可得Y在k+1時刻的狀態(tài)估計值k+1|k+1及其對應(yīng)的協(xié)方差矩陣Py,k+1|k+1.

由上述推導(dǎo)可以看出，BVCKF算法考慮了濾波過程中出現(xiàn)兩個隨機(jī)變量的情況，為同時處理目標(biāo)與外輻射源的狀態(tài)估計提供了可靠的理論工具。同時采用球面-徑向容積準(zhǔn)則來處理非線性函數(shù)的多維向量積分。這種理論方法上的改進(jìn)可以從原理上降低APCL系統(tǒng)中外輻射源狀態(tài)不確定性對目標(biāo)狀態(tài)估計造成的影響。相比于以往對單變量的SA方法[9,11]和新息協(xié)方差修正算法[11]而言，在算法精度和運(yùn)算速度上都有著顯著的提升。

3 仿真分析

為了充分驗證所提算法的有效性，本文選取目標(biāo)狀態(tài)估計的均方根誤差(RMSE)和算法的單次運(yùn)行耗時作為性能評價指標(biāo)，對比分析了DP-SA算法[11]和本文所提的DP-BVCKF算法在定位跟蹤性能上的差異。由于文獻(xiàn)[11]已經(jīng)對DP-SA算法中DP階段提取目標(biāo)量測的性能進(jìn)行了較為詳細(xì)的仿真分析，本文不再贅述。

3.1 仿真參數(shù)設(shè)置

目標(biāo)、外輻射源以及觀測站的初始狀態(tài)如表1所示。

表1 仿真場景的初始狀態(tài)

目標(biāo)和外輻射源共同的過程噪聲為

Γk=diag{(5 m)2,(5 m)2,(0.5 m/s)2,(0.5 m/s)2}。

量測噪聲的協(xié)方差表示為

Rk=diag{(0.005 rad)2,(0.005 rad)2,(50 m)2,(2 m/s)2}。

目標(biāo)等效橢球體的長短軸半徑B=1 m、C=0.5 m，外輻射源輻射功率密度E=1 W/m3. 雜波的密度參數(shù)λ=0.01 m-1·rad-1，幅值參數(shù)ζ=0.5，樣本數(shù)c=6 000，系統(tǒng)的檢測概率Pd=0.95,PFT=0.05. 探測時長M=50 s，探測間隔ΔT為1 s和在21～25 s時，目標(biāo)和外輻射源做機(jī)動轉(zhuǎn)彎運(yùn)動，轉(zhuǎn)彎率ω設(shè)定為0.1 rad/s. 其余時間目標(biāo)和外輻射源的轉(zhuǎn)彎率為0 rad/s，即做勻速直線運(yùn)動；觀測站在整個探測周期內(nèi)作ω=0.1 rad/s的勻速圓周運(yùn)動。蒙特卡洛仿真次數(shù)設(shè)定為200次?；贒P-BVCKF算法的仿真場景如圖3所示。

圖3 基于DP-BVCKF算法的仿真場景示意圖

3.2 算法跟蹤性能對比分析

為了更好地驗證分析DP-SA算法和本文所提的DP-BVCKF算法在性能上的差異，本文分別基于不同的初始協(xié)方差和不同的轉(zhuǎn)彎率給出了對比仿真的結(jié)果，如圖4和圖5所示。

圖4 不同初始協(xié)方差下算法性能對比圖

從圖4和圖5中可以看出：在目標(biāo)做勻速運(yùn)動時，兩種算法都能取得很好的收斂效果，且無論處于何種情形下，DP-BVCKF算法的收斂精度都更高。這就很好地驗證了本文的理論分析，因為在DP階段兩種算法的性能效果一致;而在狀態(tài)估計階段，BVCKF在原理上具有根本性的優(yōu)勢，既能避免SA增加的額外計算量，又能實時對外輻射源狀態(tài)進(jìn)行同步的更新估計，極大地降低了外輻射源狀態(tài)不確定性的影響。但是DP-SA算法和本文所提的DP-BVCKF算法都存在一個明顯的缺陷，即對目標(biāo)機(jī)動時的跟蹤效果不夠理想，容易出現(xiàn)濾波發(fā)散的現(xiàn)象。而且，由于DP-SA算法中以擴(kuò)展卡爾曼濾波作為核心濾波工具，相比于DP-BVCKF算法而言，濾波發(fā)散現(xiàn)象更加嚴(yán)重。

對比圖4(a)、圖4(b)、圖4(c)、圖4(d)可以看出，當(dāng)初始協(xié)方差增大后，兩種算法對目標(biāo)和外輻射源的位置估計精度受影響較小。但是DP-BVCKF算法的狀態(tài)估計效果要明顯好于DP-SA算法。這是因為DP-SA算法進(jìn)行SA后，目標(biāo)狀態(tài)和外輻射源狀態(tài)之間并不獨(dú)立，對協(xié)方差進(jìn)行濾波處理時，二者狀態(tài)的不確定性會互相影響。

由圖5可以看出，當(dāng)目標(biāo)和外輻射源的轉(zhuǎn)彎率增大后，兩種算法對目標(biāo)和外輻射源的位置估計精度都明顯下降。但是DP-BVCKF算法受轉(zhuǎn)彎率的影響相對較小，這是因為兩種算法選用的核心濾波算法在原理上存在差異。仿真結(jié)果表明在處理非線性濾波跟蹤問題時，基于球面-容積準(zhǔn)則的確定性采樣近似方法性能要明顯優(yōu)于1階泰勒展開式的線性化處理方法。

表2給出了兩種算法在不同雜波密度下的運(yùn)行耗時。由表2可以看出：雜波密度的增大使得目標(biāo)狀態(tài)提取過程的耗時增加；由于DP-SA算法在SA時增加了一定的計算量，導(dǎo)致其在狀態(tài)估計環(huán)節(jié)產(chǎn)生的耗時要高于DP-BVCKF算法。

4 結(jié)論

本文針對雜波環(huán)境下空基外輻射源系統(tǒng)的定位問題，提出了一種基于DP-BVCKF的定位跟蹤算法。該方法首先利用DP算法從含有雜波的量測數(shù)據(jù)中提取出源于目標(biāo)的量測信息，然后提出BVCKF算法來降低外輻射源狀態(tài)不確定性對目標(biāo)狀態(tài)估計的影響。通過與DP-SA算法的仿真對比分析表明，本文所提算法能有效地提取出目標(biāo)的量測信息，并極大地提高目標(biāo)狀態(tài)估計的精度，在兼顧實時性的同時保持較高的準(zhǔn)確性。但是通常情況下雜波的統(tǒng)計特性是復(fù)雜多變的，本文僅考慮了雜波統(tǒng)計特性已知的情況下對目標(biāo)數(shù)據(jù)的提取問題，當(dāng)雜波統(tǒng)計特性未知時，該方法需要進(jìn)一步改進(jìn)。而且在非線性非高斯環(huán)境下，目標(biāo)狀態(tài)估計面臨的問題將會進(jìn)一步復(fù)雜化，尤其是在應(yīng)對機(jī)動目標(biāo)跟蹤，以及多目標(biāo)跟蹤時，需要對算法做出調(diào)整，這將是后續(xù)的研究重點(diǎn)。