侯利明,連峰,王偉
(1.西安交通大學智能網絡與網絡安全教育部重點實驗室,710049,西安;2.陜西省組合與智能導航重點實驗室,710068,西安;3.中國電子科技集團公司第二十研究所導航重點實驗室,710068,西安)
近20年來,隨機有限集(RFS)理論[1-2]被廣泛地應用于多目標跟蹤領域。與傳統的多目標跟蹤算法[3-5]相比,基于RFS理論的多目標跟蹤算法不需要明確量測與目標之間復雜繁瑣的對應關系,尤其適用于多目標跟蹤問題。目前,基于RFS的多目標跟蹤濾波器主要有概率假設密度(PHD)濾波器[6]、勢概率假設密度(CPHD)濾波器[7]、多目標多伯努利(MeMBer)濾波器[1]和δ-廣義標記多伯努利(δ-GLMB)濾波器[8-9]。特別地,δ-GLMB濾波器不僅在跟蹤精度方面優于其他RFS濾波器,且能自動生成目標航跡,是首個理論可證的具有精確閉式解的多目標貝葉斯跟蹤濾波器。由于δ-GLMB濾波器的遞推過程是嚴格封閉的,因此它對目標個數及其狀態的估計精度高于其他同類濾波器,但同樣面臨計算量大的問題。該濾波器的計算復雜度理論上與目標的個數呈線性關系,與量測的個數呈3次方關系[10],但在其實現過程中仍需要花費較大的計算代價用于假設排序問題。文獻[8-9]給出了該濾波器在線性高斯模型下的高斯混合(GM)實現和非線性模型下的序貫蒙特卡洛(SMC)實現方法,其中,由于SMC-δ-GLMB濾波算法受粒子數和粒子退化等問題的影響,盡管可以通過增加粒子數以獲得較高的跟蹤精度,但大量的粒子同樣意味著需要消耗更大的計算量,難以實現δ-GLMB濾波器對目標的快速跟蹤。在實際的多目標跟蹤應用中,非線性模型則更為普遍,如雷達和聲吶等量測均為非線性[11],因此研究非線性模型下δ-GLMB濾波器的多目標跟蹤問題具有非常重要的意義。
GM-δ-GLMB濾波器在線性高斯模型下具有閉式解,并且可以將GM-δ-GLMB濾波器與擴展卡爾曼(EK)[11]和無味卡爾曼(UK)[12]相結合用于解決非線性模型下的多目標跟蹤問題。由于EK算法只是簡單地將非線性系統模型進行局部線性化處理,存在嚴重的模型描述誤差,當模型的非線性程度較高時,該濾波算法可能會導致被跟蹤目標丟失,且在對非線性函數進行Taylor級數展開時,需要求解Jacobians矩陣和Hessians矩陣(二階EKF),而這2種矩陣的計算容易導致數值不穩定,且計算量較大,更有甚者,它們根本就不存在[13]。UK算法則是通過選擇一組確定的采樣點來近似非線性函數的均值和方差,精度可達到二階,但其采樣點及權系數一般由狀態維數和憑慣例或經驗而設定的參數共同來決定,很難達到較高精度的跟蹤效果。關于EK和UK濾波算法的上述限制條件,直接阻礙了δ-GLMB濾波器向非線性多目標跟蹤應用領域的推廣。
本文將積分卡爾曼(QK)[14-15]非線性濾波算法與GM-δ-GLMB濾波器相結合,提出了具有更加嚴謹的數學理論基礎且適用性更強的QK-GM-δ-GLMB濾波器。為了驗證該算法的有效性,通過在不同的雜波密度和檢測概率條件下的仿真實驗,將所提算法的濾波精度和實時性能分別與EK-GM-δ-GLMB、UK-GM-δ-GLMB和SMC-δ-GLMB 3種濾波器做了較詳細的對比。實驗結果表明:QK-GM-δ-GLMB濾波器的跟蹤精度是4種濾波器中最優的,盡管其時間消耗有所增加,但與另外兩種δ-GLMB濾波器的GM實現算法均保持在同一量級,這在許多實際應用中是可以接受的。
假設每個目標均服從如下的非線性運動和量測模型,即
xk=f(xk-1)+wk-1
(1)
zk=h(xk)+vk
(2)
式中:k為觀測時刻;xk為k時刻目標的狀態向量;zk為k時刻傳感器的觀測向量;f(·)和h(·)分別為已知的非線性狀態轉移函數和非線性觀測函數;wk-1和vk分別表示相互獨立的過程噪聲和量測噪聲,并假設它們服從均值為零、協方差矩陣分別為Qk-1和Rk的高斯分布。
積分卡爾曼濾波算法的主要思想是基于Gauss -Hermite數值積分規則求取一組帶權重的積分點,利用這些積分點求取密度函數的均值和協方差估計[14-15]。假設向量x~N(x;0,Inx),則nx維的積分規則可由下式表示
(3)

通過構造一個具有零對角元素的對稱三對角矩陣J來求取其特征值,其中,εl表示矩陣J的第l個特征值,且該矩陣其他元素可由下式求得
(4)

根據上述所求取的一組帶權重積分點,可以推導出基于QK的δ-GLMB濾波器的高斯混合實現過程,具體預測步驟如下。假設k時刻的多目標后驗概率密度具有如下的δ-GLMB RFS形式[9]
(5)
式中:ξ為直到k時刻的所有關聯映射對應于標簽集合I的關聯歷史,ξ=(θ1,…,θk),θk表示k時刻航跡標簽到傳感器量測的關聯映射;X表示帶標簽的多伯努利狀態集合;數對(I,ξ)表示當前時刻所有航跡的標簽集合I∈F(L)與直到當前時刻的標簽-量測關聯歷史ξ∈Ξ之間的一種對應關系;w(I,ξ)表示該假設為真的概率;F(L)表示由標簽空間L的子集所構成的集合;離散空間Ξ表示航跡標簽到量測的關聯過程。直到k時刻的所有關聯歷史ξ條件下標簽為l對應目標航跡的廣義標記多伯努利分量的空間分布密度函數由J(ξ)(l)個高斯項組成
(6)
由δ-GLMB分布的共軛先驗性[8]可知,k+1時刻的多目標預測分布密度也具有δ-GLMB形式,即

(7)

(8)
式中:pB(x,l)表示新生目標的狀態空間密度;HL(l)與HB(l)分別表示存活目標和新生目標的示性函數。標簽為l∈L的存活航跡的空間分布密度具有如下的高斯混合形式
(9)

(10)
(11)
(12)
為簡單起見,假設目標航跡存活概率獨立于標簽l和狀態x,即pS(x,l)=pS,航跡消亡概率為qS(x,l)=1-pS(x,l)=1-pS。
設新生目標概率密度可表示為如下的LMB RFS形式
(13)
(14)

(15)

(16)
式中:I+?L+=L∪B表示預測目標的標簽空間,而新生目標航跡的假設wB(J)則與具體的應用場景相關。
用求得的積分點bl和相應權重wl來近似式(9)中第j個高斯項的均值和協方差,具體實現步驟如下。
步驟1對式(6)中組成標簽為l的航跡中的第j個高斯項的協方差矩陣進行Cholesky分解
(17)
步驟2計算組成標簽l的第j個高斯項的積分點
(18)
步驟3非線性狀態方程傳遞積分點
(19)
步驟4計算預測狀態均值
(20)
步驟5計算狀態預測協方差矩陣
(21)

(22)
則更新后的后驗密度仍然為δ-GLMB形式
δI+(L(X))[p(ξ,θ)(·|Z)]X}
(23)
式中Θ表示航跡標簽到量測關聯θ的映射空間,即L+→{0,1,…,|Z|},0表示漏檢。此時的后驗標簽集合與預測航跡標簽集合相等,仍為I+。更新假設的概率為

(24)
令k+1時刻標簽l的航跡關聯映射為θ(l),若θ(l)>0,表示該航跡被傳感器收到的量測集合zθ(l)更新;若θ(l)=0,則表示沒有量測與之關聯,該航跡漏檢??偟牧繙y更新廣義似然函數可表示為
(25)
假設檢測概率獨立于標簽l和狀態x,即pD(x,l)=pD,漏檢概率qD(x,l)=1-pD(x,l)=1-pD。
對于更新航跡,標簽為l的航跡關聯映射滿足θ(l)>0,式(25)的廣義似然函數可具體表示為
(26)
將式(22)與式(26)相乘并進行歸一化,即可得航跡l更新后的δ-GLMB空間分布密度

(27)
(28)
歸一化常數
(29)
(30)
式中:q(ξ,j)(z(θ(l)),l)表示目標(x,l)生成量測z(θ(l))的似然概率;κ(z(θ(l)))表示量測過程中的雜波密度函數。
對于漏檢航跡,因沒有量測與之關聯,故θ(l)=0,此時的廣義似然函數可重新表示為
(31)

(32)
將式(22)與式(31)相乘并除以式(32),可得標簽l的航跡漏檢時的空間分布密度,此時的假設更新概率、狀態均值和狀態協方差矩陣分別對應于該航跡的預測空間密度式(22)中的各參數,即
(33)
(34)
(35)

步驟1對式(22)中組成航跡l的第j個高斯項的協方差矩陣進行Cholesky分解
(36)
步驟2計算標簽l的第j個高斯項的積分點
(37)
步驟3非線性量測方程傳遞積分點
(38)
步驟4計算預測量測均值
(39)
步驟5計算量測預測協方差矩陣
S(ξ,j)=
(40)
步驟6計算互協方差矩陣
(41)
步驟7計算卡爾曼增益
K(ξ,j)=G(ξ,j)[S(ξ,j)]-1
(42)
步驟8均值更新
(43)
步驟9協方差矩陣更新
(44)
在二維平面內,設觀測區域范圍為[0°,180°]×[0 m,2 000 m],觀測時長為50 s,傳感器采樣間隔T=1 s,最多可能同時有8個目標隨機出現在該區域的任何位置,且目標出現和消失的時刻都是不確定的。目標真實運動軌跡如圖1所示。

圖1 極坐標下多目標的真實運動軌跡
目標的運動模型為CT模型[4],且觀測數據存在角度和距離的觀測噪聲,k時刻目標的狀態向量xk=[px,k,vx,k,py,k,vy,k,ωk]T由位置分量(px,k,py,k)、速度分量(vx,k,vy,k)和轉彎速率ωk組成。CT模型的狀態轉移方程可表示為
xk=F(ωk-1)xk-1+Gkwk-1
(45)

F(ωk-1)=
(46)
(47)
含方位角和徑向距離噪聲的非線性量測方程為
(48)

目標的存活概率pS,k=0.99,檢測概率pD,k=0.98。雜波模型服從強度為κk=λcVu(z)的泊松RFS,其中λc表示平均每幀雜波的個數,V為傳感器監控區域的面積,u(·)為該區域內的均勻概率分布。
采用剪切和合并技術限制高斯項數目的增長,設航跡的截斷概率閾值T′=10-5,高斯項權重的閾值T″=10-5,合并閾值U=4,航跡最大值Mmax=300,每條航跡對應高斯項個數的最大值Jmax=100,其具體剪切和合并過程與文獻[16]所述一致。
實驗通過100次蒙特卡洛(MC)仿真驗證了所提算法的有效性,采用最優子模式分配(OSPA)距離DOSPA[17]來評價幾種多目標跟蹤算法的性能。本實驗中,將OSPA表達式的誤差調節參數c和階次參數p分別設置為100、2。與已有的EK-GM-δ-GLMB、UK-GM-δ-GLMB和SMC-δ-GLMB濾波算法關于平均OSPA距離和平均單步運行時間在不同檢測概率和不同的雜波強度條件下做了比較,其中,在SMC實現過程中,每個假設航跡被分配的最大粒子數為10 000,最小粒子數為2 000。
仿真實驗運行環境軟件和硬件參數為Matlab R2017b,Windows 7 64 bit,Intel Core i5-4570 CPU 3.20 GHz,RAM 4.00 GB。
圖2給出了4種算法在檢測概率pD為0.98、雜波密度λc為3.2×10-3(rad·m)-1的實驗場景中經100次MC仿真實驗后DOSPA隨時間變化曲線的對比。從圖2可以看出,4種算法的DOSPA具有相同的變化趨勢,說明在存在目標消失和新生目標出現的場景中,它們均可以準確跟蹤目標,但在跟蹤精度上有所差別。

圖2 100次MC仿真實驗后DOSPA的對比
由圖2可知,各濾波器跟蹤精度依QK-GM-δ-GLMB、EK-GM-δ-GLMB、UK-GM-δ-GLMB、SMC-δ-GLMB順序依次降低,具體原因分析如下。
(1)在δ-GLMB濾波器的更新過程中涉及到量測信息與航跡標簽之間的關聯問題,而這種關聯問題存在較大的不確定性,隨著量測信息的不斷增加,量測-標簽之間的關聯假設呈指數增長,造成大量粒子被用來逼近這些小權重的關聯假設,使得大量的計算資源被用于粒子的重采樣過程,從而導致跟蹤精度與計算時間效率都較低。SMC-δ-GLMB算法的跟蹤精度在很大程度上取決于粒子數的大小,隨著粒子數的增加,SMC-δ-GLMB算法的實現精度也會隨之升高,但同時意味著需要更長的計算時間,不易實現快速跟蹤。
(2)UK-GM-δ-GLMB算法的濾波參數值一般是憑經驗和慣例而得到的,針對復雜和高狀態維數的多目標跟蹤問題往往難以湊效。
(3)在本次實驗中,EK-GM-δ-GLMB算法的跟蹤精度僅次于QK-GM-δ-GLMB算法,這是因為針對弱非線性問題,EK算法是首要考慮的方法,當非線性程度較高時,很容易出現濾波發散問題,從而導致被跟蹤目標的丟失。
(4)本文所提的QK-GM-δ-GLMB濾波器具有較高的跟蹤精度,因為其所有濾波參數均是基于Gauss-Hermite積分規則得到,所以得到的參數更加準確,從而可實現較高精度的目標跟蹤。
為進一步詳細比較4種算法的目標跟蹤精度和實時性,在表1和表2中分別列出了不同場景中各濾波器經100次MC仿真實驗后的DOSPA和單步運行時間的比較結果。
表1給出了當雜波密度λc為0.8×10-3(rad·m)-1時,在檢測概率pD為0.80、0.90、0.98的條件下4種濾波器的DOSPA和單步運行時間的實驗數據。由表1可以看出,在相同的雜波密度場景下,各濾波器的DOSPA和單步運行時間都隨檢測概率的增加而減小,其中,在高檢測概率場景中,QK-GM-δ-GLMB算法的跟蹤精度最高,而在較低檢測概率場景中QK-GM-δ-GLMB與EK-GM-δ-GLMB算法同時具有最高跟蹤精度,且實時性優于SMC-δ-GLMB算法。
表2給出了在檢測概率pD為0.98的條件下,當雜波密度λc分別為0.8、1.6、3.2×10-3(rad·m)-1時4種濾波器的DOSPA和平均單步運行時間的實驗數據。由表2可以看出,在同一檢測概率條件下,各濾波器的DOSPA和單步運行時間都隨著雜波密度的增強而增加,其中,在不同的雜波強度下,QK-GM-δ-GLMB算法具有最高的跟蹤精度,同樣其實時性表現處于UK-GM-δ-GLMB算法和SMC-δ-GLMB算法之間。當雜波密度較高時,SMC-δ-GLMB算法的計算代價是巨大的。

表1 不同檢測概率條件下的DOSPA和平均單步運行時間

表2 不同雜波密度下的DOSPA和平均單步運行時間
綜上所述,本文所提的QK-GM-δ-GLMB濾波算法在不同的場景中都能夠以較高精度對目標進行準確跟蹤,而其對應的時間消耗在實際應用中處于可接受的范圍內,其適用場景也得到進一步擴展。
在非線性高斯多目標運動場景下,給出了QK-GM-δ-GLMB濾波算法的遞推過程,QK-GM-δ-GLMB濾波算法具有良好的濾波性能,是一種適應性強、估計精度高的跟蹤算法,其時間開銷在實際的多目標跟蹤應用中是完全可接受的,所以QK-GM-δ-GLMB濾波算法具有一定的工程應用價值。
由于QK-GM-δ-GLMB濾波算法在非線性場景中具有較高的跟蹤精度,下一步將致力于將其應用于機動目標、擴展目標和群目標等較復雜的多目標跟蹤問題。