徐 洋，方洋旺，伍友利，張丹旭

(空軍工程大學 航空工程學院， 陜西 西安 710038)

在目標跟蹤領域中，系統(tǒng)模型誤差不可避免，產(chǎn)生的原因主要有目標運動狀態(tài)的變化及目標運動模型的“雙非”(非線性、非高斯性)問題[1]。解決上述問題的主要方法有：提高運動模型自適應能力，如自適應機動模型[2]；改進濾波算法性能，如粒子濾波器(Particle Filter, PF)[3]；增強跟蹤模型的適用范圍以解決“雙非”問題，如高斯混合模型(Gaussian Mixture Model, GMM)[4]。由于自適應機動模型人為設定因素過多且適用性還有待檢驗，而PF算法計算復雜度過高，難以應用于戰(zhàn)時平臺，為此傳統(tǒng)方法是引入GMM對“雙非”系統(tǒng)下的目標狀態(tài)進行估計。

GMM是單一高斯概率密度函數(shù)(Probability Density Function, PDF)的擴展，理論上它可以平滑地近似任意非線性函數(shù)[5]，但其本身仍存在不足：當動態(tài)方程的噪聲項及初始狀態(tài)中不止一個采用GMM表示時，算法的計算復雜度會呈指數(shù)級增長[6]。為此學者們分別從不同角度提出解決方案，這些方法大體可概括為以下兩種：基于剪枝的刪減策略與基于聚類的合并策略。

基于剪枝的刪減策略即是基于某種測量準則剔除影響程度較低的混合項。其中最為簡單的一種即為刪掉具有較低權重的高斯混合項[7]，但該方法對于具有多峰分布的PDF近似效果并不佳[8]。為此文獻[9]提出了一種基于KL(Kullback-Leibler )散度的剪枝方法。該方法在剪枝的同時依然可保留原分布峰值特征，避免了上述問題的發(fā)生。文獻[6]更是直接利用初始狀態(tài)篩選有效的噪聲混合項，從源頭降低了混合項數(shù)量。

基于聚類的合并策略即是基于某種距離標準將相似的混合項加以合并，該削減策略得到了廣泛認可。文獻[10]通過尋找混合項中兩兩之間具有最小KL散度的混合項，然后將其合并達到削減數(shù)量的目的。但該方法每合并一次都要重新搜索所有混合項，時效性較差。文獻[11]在此基礎上，又利用了聚類合并的思想，雖提高了精度，但算法復雜度仍很高。文獻[12]更是對典型合并方法進行了系統(tǒng)性分析與比較。正是以上GMM削減方法的提出，使得GMM得以應用于眾多領域[13-17]。

為了充分利用混合項自身的分布信息，簡化整體運算流程，在保證精度的同時提高算法效率，本文提出了一種復雜度適中、近似效果較好的高斯聚類-合并方法。首先基于高斯混合項的后驗分布特性及置信區(qū)間，建立相似分布特性準則，并采用擴展積分均方誤差(Extended Integral Square Error，EISE)準則求解最優(yōu)置信區(qū)間，完成對混合項的高斯聚類；然后針對各高斯簇間存在的交叉項復用問題，利用局部最近鄰(Local Nearest Neighbor, LNN)思想將其重新分配；最后在各高斯簇中利用多元合并方法，減少下次迭代混合項數(shù)量。通過仿真分析了有效置信范圍的可行域，縮窄了置信范圍的搜索空間，并通過兩個典型場景與基于數(shù)量限制的高斯混合項[7](Gaussian Mixture Reduction based on Pruning，GMRP)、Runnall[10]及基于聚類的高斯混合項削減方法[11](Gaussian Mixture Reduction based on Clustering,GMRC)進行了對比，驗證了本文算法具有較高的效費比。

1 GMM及改進方法

設線性系統(tǒng)的狀態(tài)方程及觀測方程為

(1)

式中：k≤NK為時間序列，NK為序列長度；狀態(tài)量xk∈Rmx，其中mx為x的維度，初始狀態(tài)分布為p(x0)，轉(zhuǎn)移概率p(xk|xk-1)滿足馬爾科夫過程。觀測量zk∈Rmz，其中mz為觀測量z的維度，p(zk|xk)為系統(tǒng)似然函數(shù)；k時刻狀態(tài)量及觀測量的集合分別為x0:k={x0,…,xk}及Zk={z0,…,zk}；uk及vk分別為狀態(tài)噪聲和觀測噪聲；F為狀態(tài)轉(zhuǎn)移矩陣，H為觀測矩陣。

利用GMM表示的總體概率密度為

(2)

式中：N(x;μ(i),P(i))為多維高斯分布的PDF；N為GMM中混合項數(shù)量；ΩN代表均值、方差及權重的集合；ωi代表第i個正態(tài)分布的權值。 將初始狀態(tài)、狀態(tài)噪聲及測量噪聲中的兩個或三個用GMM加以表示，會導致在狀態(tài)更新及時間更新階段，系統(tǒng)的高斯混合項呈現(xiàn)指數(shù)級增長。 針對該點不足，在綜合考慮高斯混合項權重信息和分布特性的基礎上，本文提出如下基于相似分布特性準則的聚類-合并方法以動態(tài)削減高斯混合項。

2 混合項聚類-合并算法

2.1 高斯聚類

服從高斯分布的混合項之間無法定義為完全獨立，而且聚類問題本身在統(tǒng)計上也是無法識別的[18]，因此每種聚類方法都包含著主觀的聚類概念。 本文所提的高斯聚類方法是以混合項具有較小方差[18]為前提，依據(jù)相似分布特性準則加以聚類。 若混合項中存有較大方差項，可通過設置閾值δ將其排除在聚類項之外。 其中對于聚類依據(jù)——相似分布特性準則的定義如下：若p個具有較小方差的混合項的分布中心相互處于對方的P*置信范圍內(nèi)，則認為這些混合項在分布上具有P*相似特性，從而將其歸為P*置信范圍下的同一高斯簇。 該準則的現(xiàn)實意義即是：若高斯混合項之間分布形狀足夠相似、分布中心足夠接近，則其合并后的PDF將服從單峰分布，且總體分布特性服從類高斯分布[19]。 為了搜尋可同屬一個高斯簇的混合項，定義k時刻高斯混合項的分布中心向量Π為

(3)

定義如式(4)所示判別矩陣R用以判斷混合項之間的相似性。

(4)

其中，

(5)

為了快速完成聚類，首先按行搜索判別矩陣，并對每個混合項附上隸屬簇標簽

(6)

步驟1：將交叉的高斯簇分成以下兩種

(7)

步驟2：重新分配交叉項。

1)H1:將交叉項依LNN思想重新分配。 其過程如下：

(8)

②計算交叉項與所屬高斯簇間距離

(9)

④重復步驟①～③直至交叉項分配完。

由上述可知：P*對于聚類結(jié)果影響很大。 若P*過大，則可能使得不具有相似分布特性的混合項歸并為同一簇，使得混合項發(fā)生多樣性退化，最終導致近似多峰分布時效果變差；若P*過小，則會使得聚類過細，無法明顯減少高斯混合項個數(shù)而使得實效性達不到要求。 因此合理地選擇P*對于提高算法的效費比尤為重要。 為此，本文提出基于EISE代價函數(shù)求解P*。

2.2 基于EISE函數(shù)的置信范圍搜尋算法

為了篩選合適的P*，本文受文獻[20]中所提ISE函數(shù)啟發(fā)，通過增加正則化項，在限定高斯簇個數(shù)的同時使得近似誤差逐漸降低，得式(10)所示EISE函數(shù)。

(10)

=Jhh-2Jhr+Jrr

(11)

其中，

式(11)兩個高斯PDF可表示為

(12)

將式(12)代入到式(11)中，可得

(13)

兩個高斯PDF的乘積可以簡化為

N{x;u1,P1}·N{x;u2,P2}=αN{x;u3,P3}

(14)

其中，

α=N{u1;u2,P1+P2}

將式(14)代入式(13)，又因高斯子項的積分是在整個輸入空間，故其積分結(jié)果為1，則最后只剩α。 則式(13)可表示為

(15)

至此，代價函數(shù)等于高斯簇中混合項之間的乘積加上每一個混合項與合并后項的乘積，加上合并后混合項的概率平方，再加上一個關于高斯簇數(shù)量的懲罰函數(shù)。因高斯簇已將全部混合項劃分為多個高斯簇，而高斯簇內(nèi)混合項數(shù)量也隨著P*而變化——P*越小，混合項數(shù)量越少，這對于算法效率的提高是有幫助的，但同時也使得懲罰項變大。為了更有效地尋找最優(yōu)解，可先縮小搜索域后利用MATLAB自帶函數(shù)fmincon求解。

2.3 無偏多元素合并算法

為了降低聚類復雜度，在獲得混合項后可剔除掉連續(xù)兩步權值都非常小的混合項，此處可設閾值為0.01。之后基于文獻[20]提出的合并方法，在保證系統(tǒng)狀態(tài)無偏的前提下，對各高斯簇內(nèi)的混合項進行合并

(16)

(17)

(18)

該融合方法可以保持原混合高斯項的中心矩不變[20]，而且計算方便，被大多數(shù)基于合并思想的改進方法采用[8,10-11]。

3 算法分析

3.1 算法流程

為了清晰算法步驟，下文中用GMM分別表示初始狀態(tài)、系統(tǒng)噪聲及觀測噪聲，并采用卡爾曼濾波方法對目標狀態(tài)進行估計，具體算法流程可參考算法1。

算法1 聚類-合并算法

3.2 算法復雜度

在時間復雜度上，所提算法首先需要判斷混合項之間是否具有相似分布特性，即計算判別矩陣。而判別矩陣只需得知上三角的值即可判別所有混合項間的相似關系，故第一步的時間復雜度為O(n2)；第二步因交叉項自帶標簽，故查找交叉項并不會占用時間，但在重新分配時需判斷其與各所屬簇間的最小距離，因此時間復雜度為O(n2)；第三步因采用多元素融合方法，其時間復雜度為O(lgn)；第四步為對置信范圍采用非線性優(yōu)化方法尋優(yōu)，其時間復雜度為O(m×p×n2)，其中m為算法迭代次數(shù)，p為搜索粒子數(shù)量。由此可見，用于搜索的時間復雜度受代價函數(shù)性能、粒子初始狀態(tài)及數(shù)量的影響，并且隨著狀態(tài)維度的升高，用于搜索最優(yōu)P*值的迭代次數(shù)會呈指數(shù)級增加。因此所提算法更適用于低維度的狀態(tài)。而為了限制該步驟的時間復雜度，有效合理地縮窄搜索范圍至關重要。

針對空間復雜度，因算法需要空間存儲融合后的混合項、高斯簇，故其空間復雜度為O(n)。標準算法因在迭代過程中需要大量空間來存儲中間變量，并且每次采樣混合項數(shù)量呈指數(shù)式增長，所以在空間復雜度上，改進算法具有明顯優(yōu)勢。

4 性能分析

為了驗證算法的有效性，首先利用三種典型場景來確定置信范圍的有效搜索域，從而提高算法的搜索效率。同時為了方便圖形呈現(xiàn)，選取一維情況加以說明。最后，針對二種常見的目標運動場景對運動目標的狀態(tài)進行估計，以均方根誤差(Root Mean Square Error, RMSE)作為評價跟蹤性能的依據(jù)，并以中央處理單元(Central Processing Unit,CPU)運行時間作為評價時間復雜度的標準。

4.1 置信范圍有效區(qū)間

4.1.1 四混合項呈單峰分布情況

四個混合項的參數(shù)分別設置如下：

Ω1(1)={0.3,N(x;1,12)}

Ω1(2)={0.1,N(x;1.5,1.22)}

Ω1(3)={0.4,N(x;2,0.52)}

Ω1(4)={0.2,N(x;3,0.82)}

令P*分別為1σ,2σ,3σ,4σ時，PDF近似效果如圖1所示。由圖可知，P*在1σ～3σ時算法對單峰PDF近似效果較好，而在P*=4σ時效果明顯變差。這主要是因為當P*選擇過大時，無法等效成單峰分布的混合項被歸并成一簇，從而導致混合項多樣性退化，導致近似精度下降。

圖1 近似單峰PDFFig.1 Approximated PDF to single peak

表1為不同P*下算法效率。由表可知，當P*取2σ和3σ時近似效果幾乎相同，這說明單純的增加或減小P*并不一定能夠改善近似效果，這與當前混合項的分布特性有關。并且，隨著P*的增大，最終的混合項個數(shù)將呈下降趨勢。

表1 單峰PDF下算法效率

4.1.2 四混合項呈雙峰分布情況

四個高斯混合項的參數(shù)設置如下所示：

Ω2(1)={0.3,N(x;0.5,12)}

Ω2(2)={0.1,N(x;2,1.22)}

Ω2(3)={0.4,N(x;3.5,0.52)}

Ω2(4)={0.2,N(x;4,0.82)}

圖2為置信范圍分別為1σ、2σ、3σ及4σ時對應的PDF曲線。由圖可以看出PDF的近似效果隨著置信范圍的增加而逐漸變差。當P*=3σ或P*=4σ時，近似效果不佳，尤其在峰值處。

圖2 近似雙峰PDFFig.2 Approximated PDF to double peak

該情況下算法效率如表2所示。由表可知，在近似具有雙峰分布的PDF時，算法在降低高斯項個數(shù)方面性能較好，當P*=4σ時，其返回的高斯項個數(shù)僅為初始輸入的一半，但耗時最長。

表2 雙峰下PDF算法效率

4.1.3 八混合項呈多峰分布情況

八個混合項的參數(shù)設置如下：

Ω3(1)={0.1,N(x;0.5,1)}

Ω3(2)={0.05,N(x;2,1.22)}

Ω3(3)={0.35,N(x;3.5,0.52)}

Ω3(4)={0.1,N(x;4,0.82)}

Ω3(5)={0.1,N(x;5,0.22)}

Ω3(6)={0.2,N(x;-1,0.32)}

Ω3(7)={0.04,N(x;0,42)}

Ω3(8)={0.06,N(x;1,32)}

圖3為不同置信范圍P*下對應的近似PDF。由圖可知在近似多峰分布的PDF時，改進算法的整體近似效果較為理想，但隨著P*的增大，對于起伏處的PDF近似效果逐漸變差。算法效率如表3所示。

圖3 近似多峰PDFFig.3 Approximated PDF with multiple peaks

由表3可以看出，P*的增加使得耗時也稍有增加，但同時獲得的高斯混合項的個數(shù)也在減少，這更有利于用在遞推算法中。

通過分析上述三種典型情況，可知過大的P*會使得近似精度變差，盡管可以大幅度降低輸出的混合項個數(shù)、節(jié)省GMM迭代時間，但相應會增加聚類-合并的處理時間。

表3 多峰PDF下算法效率

在通過EISE代價函數(shù)計算的最優(yōu)置信范圍下，三種典型情況的近似性能如表4所示。

表4 最優(yōu)置信范圍下近似效果

通過上面三個典型的場景，將置信范圍P*的搜索域設定在[0.5σ,3σ]以提高算法效率。

4.2 跟蹤性能

在雷達接收系統(tǒng)中，因目標的散射特性使觀測噪聲也可能呈現(xiàn)非高斯性質(zhì)，稱之為閃爍噪聲[21]。不同于高斯噪聲，閃爍噪聲尾部較長，而在中心區(qū)域則類似高斯形狀。其概率分布可表示為

ft=(1-ε)fg+εfl

(19)

式中，ft,fg及fl分別代表閃爍噪聲、高斯噪聲及拉普拉斯噪聲的分布，其表達式分別為

(20)

(21)

其中：xs表示特征變量；ε為一個小于1的小正數(shù)，其表示為閃爍頻率，下文設ε=0.01；fg的方差σ小于fl的方差η。

其中參量的設定為：系統(tǒng)狀態(tài)噪聲的方差Q=diag[10,4,2]，x1=[30,20,10]T,P1=diag[300,100,50]，x2=[20,10,5]T,P2=diag[350,200,40]，v1=[10,5,2]T,R1=diag[100,25,4]，v2=[-10,3,1]T,R2=diag[100,9,1]。

4.2.1 圓周運動場景

在弱機動場景下，目標以0.02π rad/s的角速度作勻速圓周運動，其在100 s以內(nèi)位置RMSE的對比如圖4所示。

圖4 位置RMSEFig.4 Position RMSE

圖4中GMRC、Runnall、GMRCM、GMRP四種算法的平均位置RMSE分別為8.96 m, 9.29 m, 9.32 m以及9.78 m，由此可以看出對于弱機動目標，GMRC算法的跟蹤精度最好，其次Runnall與本文提出的GMRCM算法精度相當，而GMRP算法的精度最差。在第50 s時刻目標狀態(tài)的后驗分布及混合高斯近似的后驗分布如圖5所示。

圖5 目標狀態(tài)后驗分布Fig.5 Posterior distribution of target state

由圖5可知，因多模型及非高斯噪聲的存在，估計的位置后驗概率密度函數(shù)呈現(xiàn)非高斯性質(zhì)。各算法對于圓周運動的位置估計精度較為理想，其偏差幅度并不大。其中GMRC、Runnall以及GMRCM算法的估計精度相當，略優(yōu)于GMRP算法。

置信范圍P*隨時間的變化情況如圖6所示。

圖6 最優(yōu)置信范圍的變化Fig.6 Variation of the confidence interval

由圖6可知，針對弱機動情況下，P*大部分情況處于0.5～1.5P(1，1)之間(對應圖中4～14 m)，即說明了在目標狀態(tài)不發(fā)生突變的情況下，一般會取較大置信范圍值，這主要是由于當前目標的狀態(tài)后驗分布比較均勻，具有多峰分布的情況并不多才使得P*較大。而P*會發(fā)生不斷的跳動主要是由初始搜索點選擇的隨機性以及存在多個最優(yōu)解造成的。

圖7為間隔10個采樣時刻的高斯混合項個數(shù)的變化情況。

圖7 高斯混合項個數(shù)Fig.7 Number of Gaussian components

由圖7知，改進算法所得高斯混合項個數(shù)要明顯少于另外三種算法，這主要是由于GMRP算法與Runnall算法都是達到人為設定數(shù)量后才停止合并，而GMRC算法則還需要進行聚類操作，其數(shù)量會少于Runnall算法，但其用時則會受到影響。

針對算法的時間復雜性，分別記錄了GMRC、GMRP、GMRCM、Runnall這四種算法的運行時長，分別為0.52 s, 0.44 s, 0.38 s, 0.48 s。由此可看出本文算法由于高效的搜索方法，迭代過程中的高斯混合項數(shù)量少，耗時最短。而GMRP算法由于在刪除混合項過程中用時極短，即便在迭代過程中會耗時較長，但仍取得了較優(yōu)的時間復雜度；GMRC算法由于計算最為復雜，故其用在篩選合并項過程中耗時最長，時間復雜度較高；Runnall算法因?qū)ふ液喜㈨椣鄬唵危蕰r間復雜度稍好于GMRC算法。

由上述分析可知，本文所提算法雖然在近似精度上不如GMRC算法，但能夠達到與Runnall算法相當?shù)木?，并且要明顯強于GMRP算法。同時，在時間復雜度上，所提算法在對比算法中效率是最高的。

4.2.2 階躍機動場景

該場景下目標機動較強，加速度變化情況為

圖8 位置RMSEFig.8 Position RMSE

位置RMSE的對比情況如圖8所示。針對階躍機動情況，各算法平均位置RMSE為：GMRC為9.32 m、Runnall為9.58 m、GMRCM為10.2 m、GMRP為12.17 m。由此可知，雖然所提GMRCM算法的跟蹤精度要明顯好于GMRP算法，但相較于另兩個算法，尤其在目標發(fā)生階躍運動時間段內(nèi)，跟蹤精度仍有待提高。

第50 s處目標狀態(tài)后驗分布如圖9所示。由圖可知，因目標發(fā)生突變機動導致兩個模型的濾波狀態(tài)發(fā)生差異，使得模型的狀態(tài)后驗分布呈現(xiàn)雙峰分布。從各算法近似曲線不難看出GMM可較好地近似出非高斯分布特性，但因狀態(tài)的突變導致高斯混合子濾波器的濾波精度下降，PDF逼近程度并不高。但從位置RMSE曲線可以看出，精度仍可滿足要求，并且GMRC與Runnall的近似精度較高，GMRCM次之，而GMRP最差。

圖9 目標狀態(tài)后驗分布Fig.9 Posterior distribution of target state

P*變化情況如圖10所示。由圖可以看出，當目標發(fā)生階躍機動時，GMRCM算法的最優(yōu)置信范圍會變小以適應狀態(tài)的突變，進而獲得較優(yōu)的跟蹤性能。而對于目標運動平緩階段，最優(yōu)置信范圍則會增大來適應穩(wěn)定的狀態(tài)后驗分布。

圖10 最優(yōu)置信范圍的變化Fig.10 Change of the confidence interval

采樣點處的混合項個數(shù)如圖11所示。圖中本文所提算法GMRCM的高斯混合項個數(shù)要少于其他三種算法。盡管當目標發(fā)生階躍機動時，目標的高斯項個數(shù)還是會增加，這主要是由于機動時刻目標狀態(tài)的后驗分布呈多峰形式，使得高斯簇增加進而增多了數(shù)量，但同比于其他三種算法，仍有優(yōu)勢。

圖11 高斯混合項個數(shù)Fig.11 Number of Gaussian components

針對算法的時間復雜性，GMRC、Runnall、GMRCM、GMRP四種算法的計算耗時分別為0.67 s, 0.57 s, 0.45 s, 0.62 s。由此可見針對階躍機動情況本文算法可以獲得較高的跟蹤效率。

由以上分析可知，在跟蹤階躍機動目標時，本文所提算法在跟蹤精度上較之GMRC、Runnall算法稍有不足，但仍在可接受范圍內(nèi)。并且，本文算法的效率在四種對比算法中是最優(yōu)的。

5 結(jié)論

針對GMM估計非高斯系統(tǒng)狀態(tài)時會導致高斯混合項個數(shù)呈指數(shù)級增長問題，本文基于相似分布特性準則，提出了一種混合聚類-合并方法。通過定義相似分布特性準則，基于EISE代價函數(shù)尋求最優(yōu)置信區(qū)間以對混合項進行高斯聚類，建立具有不同分布特性的高斯簇。為了解決各高斯簇中的交叉項的復用問題，引入LNN思想對交叉項重新分配。為了削減迭代過程中混合項數(shù)量，采用并行多元素融合法將高斯簇內(nèi)的混合項合并。最后通過仿真驗證了算法在保持跟蹤精度及提高運算效率方面的有效性。