井沛良， 段宇， 韓超， 郭榮化， 寧小磊， 劉瑜

(1.中國華陰兵器試驗中心， 陜西 華陰 714200； 2.空軍軍醫(yī)大學， 陜西 西安 714299； 3.海軍航空大學， 山東 煙臺 264001)

0 引言

作為評估無控炮彈、有控導彈等攻擊類武器打擊精度和雷達、光電等偵察類裝備目標定位精度的重要評估指標之一，圓概率誤差(CEP)在國內(nèi)外武器裝備實驗與評估中得到廣泛應用[1-3]。與命中概率[4-5]、命中/定位準確度[6](偏差)以及命中/定位密集度[7-8](方差)等評估指標無需假設攻擊點/定位點概率分布模型不同，現(xiàn)有CEP計算方法高度依賴于攻擊點/定位點概率分布模型。由于大部分武器裝備的攻擊點/定位點都呈現(xiàn)出單一密集中心特性，可以認為其分布服從或近似服從高斯分布。同時，高斯分布可由其1階原點矩和2階中心矩完全表述，而1階原點矩和2階中心矩可通過矩估計方法由樣本數(shù)據(jù)計算得到，故目前可查文獻中CEP計算方法大都采用了高斯分布模型。然而在實際實驗中，攻擊點/定位點時常會呈現(xiàn)出多密集中心現(xiàn)象，此時采用高斯分布模型顯然不再適合。國家軍用標準GJB 6289—2008和GB/T 4882—2001指出，利用高斯分布模型進行CEP計算前需進行正態(tài)性檢驗，只有通過正態(tài)性檢驗的實驗數(shù)據(jù)才可以使用高斯分布模型進行CEP計算。而對于未通過正態(tài)性檢驗的實驗數(shù)據(jù)，如呈現(xiàn)多密集中心的實驗數(shù)據(jù)，目前尚無計算CEP的具體標準可以依據(jù)。

針對上述問題，本文擬采用高斯混合模型(GMM)[9-10]作為攻擊點/定位點通用概率分布模型，并使用期望最大化(EM)算法[11]迭代求解模型參數(shù)，進而使用二分法計算CEP指標。采用該處理思路原因有3個方面：1)GMM作為非線性貝葉斯濾波[12-13]的關鍵實現(xiàn)模型之一，使用多個不同高斯概率密度函數(shù)(PDF)的加權和近似描述任意PDF，可在參與加權高斯PDF數(shù)量足夠多的條件下，使得近似誤差收斂于0，因此特別適合于描述多密集點特性的攻擊點/定位點概率分布；2)GMM待估計參數(shù)個數(shù)隨著參與加權高斯PDF數(shù)量的增加呈線性增加，因此即使采用最簡單的最大化似然(ML)參數(shù)估計方法[14-15]或傳統(tǒng)網(wǎng)格搜索法，計算量也十分巨大。雖然也可以嘗試使用蟻群算法[16-17]、遺傳算法[18-19]、模擬退火算法[20-21]等現(xiàn)代優(yōu)化算法，但是考慮到目前EM算法求解GMM參數(shù)在工程中應用更為廣泛，故本文采用EM算法進行GMM參數(shù)的ML估計計算；3)GMM為其參數(shù)的指數(shù)和函數(shù)形式，其在圓域上的積分往往需轉換至極坐標系下進行，因此給定GMM參數(shù)和圓域積分概率求解CEP的解析表達式十分困難，基于此，本文擬在給定最大積分概率誤差約束條件下，采用二分法求解CEP.

1 GMM和EM算法

1.1 GMM

GMM基本思想為使用一組高斯分布PDF加權和來近似描述任意PDF，即

(1)

式中：x為觀測矢量；K為高斯分量總個數(shù)，wk為第k個高斯分量的權重系數(shù)；N(x;μk,Rk)為第k個高斯分量的PDF，且

(2)

1.2 EM算法

EM算法廣泛應用于ML估計的多參數(shù)求解問題。ML估計作為經(jīng)典估計理論的重要方法之一，由于其求解思想簡單、漸近有效，在很多工程問題中得到廣泛應用。ML估計參數(shù)的典型求解方法有很多，如：1)最直接的方法是對似然函數(shù)關于各個待求參數(shù)求解導數(shù)，進而求解微分方程或方程組。但該方法往往與似然函數(shù)的函數(shù)形式有很大關系，經(jīng)常因為函數(shù)形式復雜而難以求得解析解；2)網(wǎng)格搜索法[22]求解邏輯簡單，但計算復雜度隨參數(shù)數(shù)量增多呈指數(shù)增加；3)Newton-Raphson迭代法[23-24]可能會遇到收斂問題。當然還可以使用遺傳算法、蟻群算法和模擬退火算法等現(xiàn)代優(yōu)化算法求解ML參數(shù)估計問題，但由于EM算法在求解GMM參數(shù)時，收斂速度快且可以保證收斂，本文采用EM算法對GMM的ML參數(shù)進行計算。

假設觀測數(shù)據(jù)為{x1,x2,…,xM}，數(shù)據(jù)間相互獨立，待估計參數(shù)為矢量θ，則似然函數(shù)可以寫成(3)式形式：

(3)

由于對數(shù)函數(shù)不改變函數(shù)的單調(diào)性，因此ML函數(shù)可通過最大化對數(shù)似然函數(shù)

(4)

求得。EM算法的關鍵是引入隱含變量的概念來對估計問題進行分解，隱含變量的引入應使得估計問題變得簡單易行，設隱含變量為z(m)，則(4)式可寫為

(5)

式中：Qm(z(m))為隱含變量的PDF. (5)式的推導使用了Jensen不等式[25]。由Jensen不等式性質(zhì)可知，無論m取任意值，f(xm,z(m);θ)/Qm(z(m))都為固定常數(shù)時，(5)式中等號成立。又由于

(6)

從而可得

(7)

(7)式即為EM算法的E步，而EM算法的M步則是基于(7)式所得Qm(z(m))對估計問題進行分解，進而采用一般ML求解方法計算：

(8)

從(8)式可以發(fā)現(xiàn)，E步中計算Qm(z(m))需要事先已知θ，而M步中計算θmax需要事先已知Qm(z(m)). 因此，EM算法具體操作時往往采取迭代計算方法，即事先給定參數(shù)初始值θ0，然后代入(7)式得到Qm(z(m))，再把Qm(z(m))代入(8)式得到θmax，如此交替使用(7)式和(8)式進行迭代，直至相鄰兩次迭代得到的θmax基本不變時結束迭代，并把迭代結果作為最終θ估計值。

2 非高斯分布CEP求解具體步驟

2.1 觀測數(shù)據(jù)收集和整理

觀測數(shù)據(jù)是指攻擊類/偵察類裝備在指定坐標系下的命中/定位坐標值，通?？蓪懗墒噶縳的形式。坐標系一般選擇笛卡爾直角坐標系，坐標原點選在靶標或待偵察目標中心點。

2.2 GMM參數(shù)選擇

GMM選用(1)式所示形式。確定高斯分量總個數(shù)K是采用GMM的關鍵。理論上，可以使用分量盡可能多的GMM. 在求解出模型參數(shù)后，僅選取權重較大的高斯分量來重構PDF，舍棄那些權重較小的，這種方法尤其適合人不在回路時大數(shù)據(jù)的自動處理。對于武器裝備CEP評估，由于人往往都在回路中，高斯分量個數(shù)可以選擇目視密集中心個數(shù)。

2.3 EM算法實施步驟

使用EM算法求解GMM參數(shù)的具體計算實施步驟如下。

步驟1給定參數(shù)初始值。

取初始權重參數(shù)為

wk=1/K.

(9)

步驟2約定隱含變量及內(nèi)涵。

對于GMM來講，由模型生成觀測數(shù)據(jù)的物理內(nèi)涵可以理解為：依據(jù)權重參數(shù)wk來隨機選取高斯分量，然后依據(jù)所選取分量的PDF(即N(x;μk,Rk))來生成具體觀測。因此，隱含變量z(m)可以構建為觀測矢量xm與高斯分量的隸屬關系，即xm到底由哪個高斯分量生成這一隨機事件，而Qm(z(m))則表述了xm隸屬于第z(m)個高斯分量的概率。

步驟3依據(jù)當前參數(shù)模型，求解出Qm(z(m)).

(10)

步驟4依據(jù)E步所得Qm(z(m))，對GMM參數(shù)進行更新。

(11)

(12)

(13)

步驟5重復步驟3和步驟4，直到相鄰兩次所得GMM參數(shù)基本不變。

2.4 針對GMM的二分法CEP計算方法

(14)

式中：|·|為矢量2范數(shù)操作。顯然，直接求解出CEP是困難的，因為給定f(x)在|x|≤r約束區(qū)域的積分概率，很難寫出以r為自變量的解析表達式。然而，考慮到給定的r取任意值時，f(x)在|x|≤r約束區(qū)域的積分概率容易求解出，且p(r)為r的單調(diào)不減函數(shù)的現(xiàn)實，可以在給定p(r)精度ξ(后續(xù)仿真實驗中ξ取為1×10-6)要求的前提下，使用二分法迭代求解r，其基本步驟如下。

步驟1初始化。

步驟2二分法求解結果測試。

步驟3二分法收縮可行解區(qū)間。

[an+1,bn+1]=[(an+bn)/2,bn]；

(15)

否則取

[an+1,bn+1]=[an,(an+bn)/2]；

(16)

并令n=n+1.

步驟4如果n=N，則取CEP=(an+bn)/2，并結束迭代；否則跳轉至步驟2.

3 仿真結果與分析

3.1 算法單次實施典型仿真結果與分析

以一個簡單的二分量GMM為實施例進行仿真。模型真實參數(shù)設置如下：K=2，w1=2/5，w2=3/5，μ1=[2,2]T，μ2=[4,5]T，R1=diag([1,1]T)，R2=diag([1,2]T)；EM算法初始GMM參數(shù)設置為：K=3，w1=1/3，w2=1/3，w3=1/3，μ1=[1,1]T，μ2=[5,6]T，μ3=[7/2,7/2]T，R1=diag([3,4]T)，R2=diag([2,3]T)，R3=diag([3,3]T). 可得真實PDF曲面和EM算法初始GMM的PDF曲面分別如圖1和圖2所示。

圖1 真實PDF曲面Fig.1 Surface of true PDF

圖2 EM算法初始設置GMM的 PDF曲面Fig.2 Surface of initial GMM PDF derived by EM algorithm

依據(jù)真實PDF生成200個觀測數(shù)據(jù)，如圖3所示。針對生成觀測數(shù)據(jù)，EM算法最終迭代求解出GMM的PDF曲面如圖4所示。

圖3 依據(jù)真實PDF生成觀測Fig.3 Observations generated by true PDF

圖4 EM算法最終迭代求解出GMM的 PDF曲面Fig.4 Surface of GMM PDF derived by EM algorithm

圖4所表征的GMM中高斯分量的權重、均值、協(xié)方差陣依次為：w1=0.341 8，w2=0.351 3，w3=0.304 9，μ1=[1.867 3,1.675 5]T，μ2=[4.189 5,5.711 2]T，μ3=[3.399 7,3.436 8]T，R1=diag([1.018 8,0.881 0]T)，R2=diag([0.935 1,1.903 7]T)，R3=diag([1.144 0,1.640 5]T). 把圖4、圖2分別與圖1進行對比可以發(fā)現(xiàn)，EM算法迭代求解GMM的PDF曲面較初始設置GMM的 PDF曲面明顯更接近于真實PDF曲面。另外，將EM算法迭代前后GMM中高斯分量的均值和協(xié)方差陣分別與真實GMM中高斯分量的均值和協(xié)方差陣進行對比可以發(fā)現(xiàn)，EM算法迭代使得GMM前兩個高斯分量更接近于真實GMM的兩個高斯分量，這表明了EM算法的良好性能。

為了定量對比本文算法和傳統(tǒng)單高斯模型的性能優(yōu)劣，采用經(jīng)典Kullback-Leibler(KL)區(qū)分度[26]指標以及PDF差值絕對值積分指標對算法性能，即算法求解PDF與真實PDF的相似度，進行評價。KL區(qū)分度是衡量兩個PDF差異大小的經(jīng)典指標，KL區(qū)分度越小，表明兩個PDF越趨于相同。為了更好地展示不同算法性能細節(jié)，本文后續(xù)部分采用了KL區(qū)分度的對數(shù)值來表征不同函數(shù)之間的接近程度。對于給定的兩個PDF：f1(x)和f2(x)，f2(x)相對于f1(x)的KL區(qū)分度定義如下：

(17)

由于lg(·)為上凸函數(shù)，由Jensen不等式可知，dKL(f1(x)‖f2(x))≥0總是成立，且等號成立時有f1(x)=f2(x). 與KL區(qū)分度類似，PDF差值絕對值積分指標也可以對兩個PDF接近程度進行評價，PDF的f2(x)相對于f1(x)的差值絕對值積分定義如下：

(18)

顯然，兩個PDF越相似，其差值絕對值積分越小。當差值絕對值積分為0時，兩個PDF完全一致。

定義真實PDF為fT，隱含變量已知時PDF為fH，傳統(tǒng)方法求解PDF為fRSGM，初始GMM的PDF為fIGMM，EM算法迭代求解GMM的PDF為fRGMM，可得fRSGM、fIGMM、fRGMM分別與fT、fH這兩個參考基準之間的KL區(qū)分度如圖5所示，差值絕對值積分如圖6所示。

圖5 不同PDF之間對數(shù)KL區(qū)分度Fig.5 Logarithmic KL divergence between different PDFs

圖6 不同PDF之間差值絕對值積分Fig.6 Integral of absolute value of difference between different PDFs

注意：fRSGM、fIGMM、fT、fH并不隨著EM算法的迭代而發(fā)生變化，因此fRSGM、fIGMM與fT、fH之間的指標值僅在圖5、圖6中的第0次迭代位置處給出。而fRSGM與fT、fH之間的指標值為EM算法迭代次數(shù)的函數(shù)，因而以曲線形式展示。從圖5和圖6中可以看出，采用EM算法求解GMM參數(shù)，能夠快速、穩(wěn)定地收斂到一個顯著優(yōu)于單高斯模型參數(shù)的結果值。由于EM算法需要克服隱含參數(shù)來求解GMM參數(shù)，因此相對于fT，EM算法迭代結果更趨近于fH. 這恰恰表明了EM算法的優(yōu)異性能，當觀測矢量個數(shù)足夠多時，fT與fH足夠相似，此時EM算法便通過逼近fH來逼近fT.

如圖7所示給出了EM算法迭代過程中，初始GMM分量中心位置的收斂變化軌跡。同時，依據(jù)求解PDF計算出了不同PDF的CEP圓域，并以圓圈的形式畫出。

圖7 EM算法迭代收斂軌跡及不同算法CEP求解 結果對比示意圖Fig.7 Iterative convergence trace of EM algorithm and CEP results of different algorithms

從圖7中可以看出，fT與fH求解出的CEP基本一致，而本文方法求解CEP要比傳統(tǒng)方法求解CEP更接近真實PDF所對應的CEP.

3.2 算法多次Monte Carlo仿真性能曲面對比與分析

為了更充分對比算法性能，固定3.1節(jié)中初始GMM參數(shù)設置，并在變化真實GMM中高斯分量中心位置的同時，保持真實GMM中高斯分量協(xié)方差陣和高斯分量系數(shù)不變。真實GMM中高斯分量中心位置僅作平移變化，對每一平移變化后的真實GMM，采用Monte Carlo方法生成100組觀測數(shù)據(jù)。針對每一組數(shù)據(jù)分別使用傳統(tǒng)方法和本文方法求解出對應PDF及CEP，并計算所得PDF與真實PDF相似程度以及所得CEP與真實CEP的接近程度。對所有對比結果進行統(tǒng)計，即可得到傳統(tǒng)方法和本文方法平均性能。當GMM中高斯分量中心位置在一定區(qū)域多次進行平移時，便可得到一個平均性能的曲面。這里采用的性能指標有PDF之間的KL區(qū)分度(如圖8、圖9所示)、差值絕對值積分(如圖10、圖11所示)以及CEP的均方誤差(如圖12、圖13所示)。曲面對應的x軸坐標值和y軸坐標值分別表示GMM中高斯分量中心位置在x軸方向和y軸方向上的平移幅度。

圖8 傳統(tǒng)方法求解PDF與真實PDF之間對數(shù)KL區(qū) 分度均值曲面Fig.8 Surface of average logarithmic KL divergence between PDF derived by traditional method and true PDF

圖9 本文方法求解PDF與真實PDF之間對數(shù)KL區(qū) 分度均值曲面Fig.9 Surface of average logarithmic KL divergence between PDF derived by the proposed method and the true PDF

圖10 傳統(tǒng)方法求解PDF與真實PDF之間差值 絕對值積分的均值曲面Fig.10 Integral mean value of absolute value of difference between PDF derived by traditional method and true PDF

圖11 本文方法求解PDF與真實PDF之間差值絕對值 積分的均值曲面Fig.11 Integral mean value of absolute value of difference between PDF derived by the proposed method and true PDF

圖12 傳統(tǒng)方法求解CEP均方誤差曲面Fig.12 Mean square error of CEP derived by the traditional method

圖13 本文方法求解CEP均方誤差曲面Fig.13 Mean square error of CEP derived by the proposed method

從圖8、圖9、圖10和圖11可以看出，本文方法在變化參數(shù)的總體范圍內(nèi)，相較于傳統(tǒng)方法能夠取得更接近于真實PDF的PDF估計。由圖12和圖13可知，圖12曲面積分值為68.376 3，圖13曲面積分值為5.311 6，二者對比結果表明，本文方法求解CEP的均方誤差顯著小于傳統(tǒng)方法求解CEP的均方誤差。這說明，本文方法相對于傳統(tǒng)方法具有十分優(yōu)異的CEP估計性能，十分適用于觀測數(shù)據(jù)非高斯分布時的CEP估計問題。

4 結論

本文分析了傳統(tǒng)CEP求解方法在處理觀測數(shù)據(jù)呈現(xiàn)非高斯分布時存在的缺陷，并基于GMM和EM算法，提出了觀測數(shù)據(jù)服從任意分布時的CEP求解新方法。仿真結果表明，無論從KL區(qū)分度，還是從PDF差值絕對值積分來分析，所提方法求解出的PDF相較于傳統(tǒng)方法明顯能夠更接近真實PDF. 最終求解CEP的均方誤差結果也揭示了本文方法的優(yōu)異性能。