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

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

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

0 引言

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

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

1 GMM和EM算法

1.1 GMM

GMM基本思想為使用一組高斯分布PDF加權(quán)和來(lái)近似描述任意PDF，即

(1)

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

(2)

1.2 EM算法

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

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

(3)

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

(4)

求得。EM算法的關(guān)鍵是引入隱含變量的概念來(lái)對(duì)估計(jì)問(wèn)題進(jìn)行分解，隱含變量的引入應(yīng)使得估計(jì)問(wèn)題變得簡(jiǎn)單易行，設(shè)隱含變量為z(m)，則(4)式可寫為

(5)

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

(6)

從而可得

(7)

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

(8)

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

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

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

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

2.2 GMM參數(shù)選擇

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

2.3 EM算法實(shí)施步驟

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

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

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

wk=1/K.

(9)

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

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

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

(10)

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

(11)

(12)

(13)

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

2.4 針對(duì)GMM的二分法CEP計(jì)算方法

(14)

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

步驟1初始化。

步驟2二分法求解結(jié)果測(cè)試。

步驟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，并結(jié)束迭代；否則跳轉(zhuǎn)至步驟2.

3 仿真結(jié)果與分析

3.1 算法單次實(shí)施典型仿真結(jié)果與分析

以一個(gè)簡(jiǎn)單的二分量GMM為實(shí)施例進(jìn)行仿真。模型真實(shí)參數(shù)設(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ù)設(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). 可得真實(shí)PDF曲面和EM算法初始GMM的PDF曲面分別如圖1和圖2所示。

圖1 真實(shí)PDF曲面Fig.1 Surface of true PDF

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

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

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

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

圖4所表征的GMM中高斯分量的權(quán)重、均值、協(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進(jìn)行對(duì)比可以發(fā)現(xiàn)，EM算法迭代求解GMM的PDF曲面較初始設(shè)置GMM的 PDF曲面明顯更接近于真實(shí)PDF曲面。另外，將EM算法迭代前后GMM中高斯分量的均值和協(xié)方差陣分別與真實(shí)GMM中高斯分量的均值和協(xié)方差陣進(jìn)行對(duì)比可以發(fā)現(xiàn)，EM算法迭代使得GMM前兩個(gè)高斯分量更接近于真實(shí)GMM的兩個(gè)高斯分量，這表明了EM算法的良好性能。

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

(17)

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

(18)

顯然，兩個(gè)PDF越相似，其差值絕對(duì)值積分越小。當(dāng)差值絕對(duì)值積分為0時(shí)，兩個(gè)PDF完全一致。

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

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

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

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

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

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

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

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

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

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

圖9 本文方法求解PDF與真實(shí)PDF之間對(duì)數(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與真實(shí)PDF之間差值 絕對(duì)值積分的均值曲面Fig.10 Integral mean value of absolute value of difference between PDF derived by traditional method and true PDF

圖11 本文方法求解PDF與真實(shí)PDF之間差值絕對(duì)值 積分的均值曲面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)方法能夠取得更接近于真實(shí)PDF的PDF估計(jì)。由圖12和圖13可知，圖12曲面積分值為68.376 3，圖13曲面積分值為5.311 6，二者對(duì)比結(jié)果表明，本文方法求解CEP的均方誤差顯著小于傳統(tǒng)方法求解CEP的均方誤差。這說(shuō)明，本文方法相對(duì)于傳統(tǒng)方法具有十分優(yōu)異的CEP估計(jì)性能，十分適用于觀測(cè)數(shù)據(jù)非高斯分布時(shí)的CEP估計(jì)問(wèn)題。

4 結(jié)論

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