EM算法在閃爍噪聲分布參數(shù)估計(jì)中的應(yīng)用

李廣柱 陳威兵 李 瑋 張剛林 蘇 鋼

(長沙學(xué)院 長沙 410022)

0 引言

近距離跟蹤雷達(dá)，目標(biāo)回波具有角閃爍特性，如果處理不當(dāng)，會(huì)增加雷達(dá)跟蹤目標(biāo)的誤差，甚至?xí)G失目標(biāo)。目標(biāo)角閃爍噪聲的分布是非高斯的，具有較長的拖尾[1-5]。在實(shí)際應(yīng)用時(shí)，常采用高斯分布和較小概率出現(xiàn)的大方差拉普拉斯分布噪聲混合進(jìn)行建模，該分布模型是研究跟蹤濾波算法的重要基礎(chǔ)[6-7]。

本文專門研究角閃爍噪聲的分布參數(shù)估計(jì)，基于混合分布模型，采用EM算法估計(jì)分布參數(shù)，給出了分布參數(shù)估計(jì)的理論推導(dǎo)，列出了參數(shù)估計(jì)的計(jì)算流程，并進(jìn)行了仿真驗(yàn)證。理論推導(dǎo)和仿真結(jié)果表明，采用EM算法可以估計(jì)得到混合分布模型的參數(shù)。

1 閃爍噪聲分布參數(shù)的EM估計(jì)算法

1.1 角閃爍噪聲的混合分布模型

角閃爍噪聲，一般認(rèn)為服從混合分布模型為

f(x;σ,η,ε)=(1-ε)fg(x)+εfl(x)

(1)

式(1)中fg(x)、fl(x)分別表示高斯分布和拉普拉斯分布的概率分布函數(shù)，ε?1是觀測(cè)噪聲服從拉普拉斯分布的概率。其中

(2)

(3)

式(2)、式(3)中σ、η分別為高斯分布、拉普拉斯分布的參數(shù)。

1.2 EM分布參數(shù)估計(jì)算法

假設(shè)有采樣序列x=[x0,x1,…,xN-1]T，服從上述混合分布模型，可知序列x的概率分布為

(4)

式(4)中的α為服從[0,1]分布的隨機(jī)數(shù)，且α=1的概率為ε。可知似然比函數(shù)為

L(θ)=lnf(x|σ,η,ε)

(5)

似然比函數(shù)的數(shù)學(xué)期望為

(6)

式(6)中的Θ為待估計(jì)參數(shù)的集合，Θ(k)表示第k次迭代時(shí)的參數(shù)值；p(xn|g;θ)、p(xn|l;θ)分別表示在參數(shù)θ=[σ,η,ε]T的條件下，xn分別服從高斯分布和拉普拉斯分布的概率；Pg、Pl分別表示樣本服從高斯分布和拉普拉斯分布的概率；P(g|xn)和P(l|xn)表示當(dāng)采樣值為xn的時(shí)候，分別服從高斯分布和拉普拉斯分布的概率。將式(1)帶入式(6)可得混合模型條件下Q(Θ;Θ(k))的表達(dá)式?？芍植紖?shù)θ=[σ,η,ε]T的似然估計(jì)在Q(Θ;Θ(k))的極值處取到。求Q(Θ;Θ(k))對(duì)θ=[σ,η,ε]T的導(dǎo)數(shù)可得

(7)

令式(7)等于0，可得θ=[σ,η,ε]T的似然估計(jì)為

(8)

為了完成迭代的運(yùn)算，還需引入P(g|xn)和P(l|xn)的更新

(9)

利用式(8)和式(9)，可以通過迭代的方式得到θ=[σ,η,ε]T的估計(jì)，下面介紹具體的實(shí)現(xiàn)算法。

2 分布參數(shù)估計(jì)的算法流程

基于式(8)和式(9)，可得到采用EM算法估計(jì)閃爍噪聲的分布參數(shù)的計(jì)算流程如下：

1)Step1，設(shè)定k=1時(shí)，分布參數(shù)的初始值θ(1)=[σ(1),η(1),ε(1)]T，可容忍的迭代次數(shù)K和可接受的迭代誤差δ；

2)Step2，對(duì)于第k+1次迭代，計(jì)算：

(10)

3)Step3，按式(9)計(jì)算P(g|xn)(k+1)和P(l|xn)(k+1)；

4)Step4，按式(8)計(jì)算θ(k+1)=[σ(k+1),η(k+1),ε(k+1)]T；

3 仿真實(shí)驗(yàn)

在實(shí)際應(yīng)用的時(shí)候，可以采集近距離跟蹤雷達(dá)在跟蹤目標(biāo)過程中的觀測(cè)噪聲，采用EM算法對(duì)角閃爍噪聲的分布參數(shù)進(jìn)行估計(jì)。本文采用仿真的方法驗(yàn)證EM算法的有效性，改變仿真時(shí)采用的角閃爍噪聲的分布參數(shù)，通過對(duì)比仿真用到的分布參數(shù)和EM算法的估計(jì)結(jié)果，可以驗(yàn)證算法的有效性。

選取參數(shù)σ2=1、η=3、ε=0.05，仿真得到長度為4096點(diǎn)的角閃爍噪聲序列。設(shè)置分布參數(shù)的初始值為：σ2(1)=0.5、η(1)=1、ε(1)=0.1，可容忍的迭代次數(shù)K=500；可接受的迭代誤差δ=10-7。采用EM迭代算法進(jìn)行估計(jì)，得到的估計(jì)結(jié)果如圖 1所示。圖 1給出的是相對(duì)誤差，故縱坐標(biāo)的單位是%。

圖1 EM算法估計(jì)噪聲分布參數(shù)的相對(duì)誤差

由圖 1可見，迭代40次左右即可達(dá)到較好的效果。為了驗(yàn)證EM算法估計(jì)的可重復(fù)性，本文模擬了500組相同分布參數(shù)的角閃爍噪聲序列，利用EM算法估計(jì)得到分布參數(shù)：σ2、η、ε的頻數(shù)分布直方圖，如圖 2到圖 4所示。

圖2 參數(shù)σ2估計(jì)值的頻數(shù)分布直方圖

圖3 參數(shù)η估計(jì)值的頻數(shù)分布直方圖

圖4 參數(shù)ε估計(jì)值的頻數(shù)分布直方圖

通過圖2到圖4可以發(fā)現(xiàn)，η，特別是ε的估計(jì)值有較長的拖尾，為此對(duì)σ2、η和ε的估計(jì)值進(jìn)行Jarque-Bera分布檢驗(yàn)，結(jié)果如表 1所示。

表1 估計(jì)值的Jarque-Bera檢驗(yàn)

表1采用的置信度為0.05，表中h表示是否接受樣本服從高斯分布：h=0，表示樣本服從高斯分布；h=1，表示樣本不服從高斯分布，由表 1可以看出，σ2的估計(jì)值服從高斯分布，η和ε不服從高斯分布。表中p為Jarque-Bera檢驗(yàn)統(tǒng)計(jì)量，p小于置信度，說明估計(jì)值不服從高斯分布，由表 1可見，η、ε的Jarque-Bera檢驗(yàn)統(tǒng)計(jì)量皆小于置信度，其中ε的Jarque-Bera檢驗(yàn)統(tǒng)計(jì)量很小，這與圖 4是吻合的。

表2給出σ2、η和ε估計(jì)值的前兩階矩的計(jì)算結(jié)果。

表2 500次估計(jì)值的均值和標(biāo)準(zhǔn)差

通過表 2可以發(fā)現(xiàn)，采用EM算法估計(jì)混合分布模型的參數(shù)σ2和ε的效果較好，參數(shù)η的估計(jì)標(biāo)準(zhǔn)差略大，這是由于ε較小，模擬得到的角閃爍噪聲序列中服從拉普拉斯分布的噪聲點(diǎn)少，樣本數(shù)較少造成的。因此，可以通過提高總的觀測(cè)樣本來提高參數(shù)η的估計(jì)精度。

4 結(jié)束語

角閃爍噪聲是影響近距離跟蹤雷達(dá)性能的重要因素，掌握它的概率分布是改善目標(biāo)跟蹤性能的重要依據(jù)。本文提出采用EM算法估計(jì)角閃爍噪聲的分布參數(shù)，推導(dǎo)了理論模型，列出了計(jì)算流程，并進(jìn)行了仿真。通過仿真發(fā)現(xiàn)，EM算法能夠估計(jì)得到角閃爍噪聲的分布參數(shù)，其中σ2的估計(jì)值服從高斯分布，η和ε不服從高斯分布；計(jì)算結(jié)果表明，σ2、η和ε估計(jì)值的均值與真值偏差小，綜合分析可以認(rèn)為算法效果良好。下一步可以通過半實(shí)物仿真和實(shí)物仿真，估計(jì)得到角閃爍噪聲的分布參數(shù)。