基于高斯混合模型和期望最大化算法的非高斯分布圓概率誤差估計方法研究

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

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

0 引言

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

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

1 GMM和EM算法

1.1 GMM

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

(1)

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

(2)

1.2 EM算法

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

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

(3)

由于對數函數不改變函數的單調性，因此ML函數可通過最大化對數似然函數

(4)

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

(5)

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

(6)

從而可得

(7)

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

(8)

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

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

2.1 觀測數據收集和整理

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

2.2 GMM參數選擇

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

2.3 EM算法實施步驟

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

步驟1給定參數初始值。

取初始權重參數為

wk=1/K.

(9)

步驟2約定隱含變量及內涵。

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

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

(10)

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

(11)

(12)

(13)

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

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

(14)

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

步驟1初始化。

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

步驟3二分法收縮可行解區間。

[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為實施例進行仿真。模型真實參數設置如下：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參數設置為：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

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

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

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

圖4所表征的GMM中高斯分量的權重、均值、協方差陣依次為：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進行對比可以發現，EM算法迭代求解GMM的PDF曲面較初始設置GMM的 PDF曲面明顯更接近于真實PDF曲面。另外，將EM算法迭代前后GMM中高斯分量的均值和協方差陣分別與真實GMM中高斯分量的均值和協方差陣進行對比可以發現，EM算法迭代使得GMM前兩個高斯分量更接近于真實GMM的兩個高斯分量，這表明了EM算法的良好性能。

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

(17)

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

(18)

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

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

圖5 不同PDF之間對數KL區分度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算法的迭代而發生變化，因此fRSGM、fIGMM與fT、fH之間的指標值僅在圖5、圖6中的第0次迭代位置處給出。而fRSGM與fT、fH之間的指標值為EM算法迭代次數的函數，因而以曲線形式展示。從圖5和圖6中可以看出，采用EM算法求解GMM參數，能夠快速、穩定地收斂到一個顯著優于單高斯模型參數的結果值。由于EM算法需要克服隱含參數來求解GMM參數，因此相對于fT，EM算法迭代結果更趨近于fH. 這恰恰表明了EM算法的優異性能，當觀測矢量個數足夠多時，fT與fH足夠相似，此時EM算法便通過逼近fH來逼近fT.

如圖7所示給出了EM算法迭代過程中，初始GMM分量中心位置的收斂變化軌跡。同時，依據求解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要比傳統方法求解CEP更接近真實PDF所對應的CEP.

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

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

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

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

圖10 傳統方法求解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 傳統方法求解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可以看出，本文方法在變化參數的總體范圍內，相較于傳統方法能夠取得更接近于真實PDF的PDF估計。由圖12和圖13可知，圖12曲面積分值為68.376 3，圖13曲面積分值為5.311 6，二者對比結果表明，本文方法求解CEP的均方誤差顯著小于傳統方法求解CEP的均方誤差。這說明，本文方法相對于傳統方法具有十分優異的CEP估計性能，十分適用于觀測數據非高斯分布時的CEP估計問題。

4 結論

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