非高斯碼元檢測(cè)的馬爾可夫鏈蒙特卡洛算法*

馮士民，周穗華，應(yīng)文威

非高斯碼元檢測(cè)的馬爾可夫鏈蒙特卡洛算法*

馮士民1，周穗華1，應(yīng)文威2

(1. 海軍工程大學(xué) 兵器工程系, 湖北 武漢 430033； 2. 中國(guó)人民解放軍91635部隊(duì), 北京 102249)

針對(duì)實(shí)際超低頻接收機(jī)不僅受非高斯噪聲的影響，還受接收機(jī)內(nèi)部和外部環(huán)境中高斯噪聲影響的問(wèn)題，對(duì)噪聲采用非高斯分布和高斯分布的混合模型建模，根據(jù)混合模型的性質(zhì)，設(shè)計(jì)了一種利用馬爾可夫鏈蒙特卡洛方法的超低頻信號(hào)碼元盲檢測(cè)算法。盲檢測(cè)算法在貝葉斯層次模型下，采用Gibbs抽樣和M-H抽樣更新參數(shù)，同步估計(jì)信道衰落系數(shù)和噪聲模型參數(shù)，并實(shí)現(xiàn)對(duì)信號(hào)碼元的檢測(cè)。算法迭代效率快、精度高。通過(guò)與最優(yōu)檢測(cè)算法性能比較，盲檢測(cè)算法性能優(yōu)異，對(duì)超低頻信號(hào)接收具有重要的現(xiàn)實(shí)意義。

非高斯噪聲；盲檢測(cè)；高斯尺度混合；混合模型

(1.DepartmentofWeaponryEngineering,NavalUniversityofEngineering,Wuhan430033,China;

2.ThePLAUnitof91635,Beijing102249,China)

超低頻通信系統(tǒng)中，受雷電等產(chǎn)生的大氣噪聲的影響，噪聲往往具有明顯的非高斯特性。傳統(tǒng)的線性接收機(jī)或高斯接收機(jī)在非高斯噪聲環(huán)境下性能急劇惡化，嚴(yán)重影響通信系統(tǒng)的正常工作，為了實(shí)現(xiàn)超低頻信號(hào)碼元的最佳檢測(cè)，需要對(duì)噪聲進(jìn)行建模和參數(shù)估計(jì)。ConteE，BuzziS和LopsM等研究了采用高斯尺度混合(GaussianScaleMixture，GSM)分布模型對(duì)噪聲建模[1-3]，并實(shí)現(xiàn)信號(hào)碼元的檢測(cè)。但是，實(shí)際中的接收機(jī)，不僅受到非高斯噪聲的影響，同時(shí)也會(huì)受到接收機(jī)內(nèi)部和外部環(huán)境中高斯噪聲的影響[4-7]。所以本文在對(duì)噪聲建模時(shí)，采用了GSM分布和高斯分布的混合模型，這對(duì)于實(shí)際中實(shí)現(xiàn)對(duì)信號(hào)碼元的最優(yōu)檢測(cè)有更好的適用性。同時(shí)，在信號(hào)碼元檢測(cè)方法上，采用傳統(tǒng)的最大期望算法或最大似然估計(jì)設(shè)計(jì)檢測(cè)算法[8-10]。由于該模型參數(shù)較多，采用上述方法時(shí)，概率密度計(jì)算難度大，精度也不夠理想。馬爾可夫鏈蒙特卡洛(MarkovChainMonteCarlo，MCMC)算法，能夠解決具有高維度且形式復(fù)雜的未知參數(shù)的后驗(yàn)概率計(jì)算問(wèn)題，是一種在統(tǒng)計(jì)計(jì)算中性能優(yōu)越的方法[11-12]。通過(guò)設(shè)計(jì)MCMC方法來(lái)實(shí)現(xiàn)信號(hào)碼元在混合噪聲模型下的盲檢測(cè)算法，迭代收斂快，精度高，具有明顯的優(yōu)勢(shì)。

1 噪聲模型

1.1 GSM模型

GSM模型是用一個(gè)高斯分布變量和一個(gè)隱含的尺度隨機(jī)變量相乘建立的模型，如式(1)所示。

(1)

n～N(0,γδ2), γ～I(xiàn)G(μ,λ)

(2)

其中，δ2為高斯分布的方差，μ和λ分別為逆高斯分布的均值和形狀參數(shù)。通過(guò)計(jì)算隨機(jī)變量n的特征函數(shù)表達(dá)式，得出：

(3)

可以看出，對(duì)于同一個(gè)分布n，λ，μ和δ2滿足特定比例關(guān)系，這表明，μ可以歸一化到高斯分布方差δ2中。為了更直觀地表示，固定參數(shù)μ為常數(shù)，這與固定逆高斯分布的平均功率為1是一致的。

1.2 噪聲模型

為了對(duì)噪聲進(jìn)行準(zhǔn)確估計(jì)并進(jìn)行信號(hào)碼元檢測(cè)，采用高斯尺度混合分布和高斯分布的混合模型作為噪聲模型：

p(n)=ω1f(n;δ2,μ,λ)+ω2N(n;0,σ2)

(4)

其中，ω1和ω2為權(quán)重因子，滿足ω1+ω2=1。

根據(jù)1.1節(jié)所述性質(zhì)，可以固定μ=5，則混合模型需要估計(jì)的參數(shù)為{δ2,λ,σ2,ω}，ω={ω1,ω2}為二維向量。

對(duì)于噪聲數(shù)據(jù)集N={ni,i=0,1,…,N}，引入標(biāo)簽變量Z={zi,i=0,1,…,N}對(duì)ni進(jìn)行分組，其中zi的取值為：當(dāng)ni屬于高斯尺度混合分布時(shí)，zi的取值為1；當(dāng)ni屬于高斯分布時(shí)，zi的取值為2。zi為獨(dú)立隨機(jī)變量，并且滿足

p(zi=j)=ωj(j=1,2)

(5)

根據(jù)式(3)、式(4)和式(5)，可以得出噪聲數(shù)據(jù)集的等效表達(dá)式：

(6)

2 信號(hào)碼元檢測(cè)模型

2.1 接收信號(hào)模型

建立平坦衰落信道下的超低頻信號(hào)碼元接收模型，基于以下3點(diǎn)假設(shè)：①假設(shè)信號(hào)碼元在超低頻頻段上進(jìn)行傳輸；②信道中的噪聲為加性、混合模型噪聲；③信道是平坦衰落的，并假設(shè)在處理時(shí)間窗內(nèi)衰落系數(shù)不變。

設(shè)S1,…,SN為發(fā)射的信號(hào)碼元，每個(gè)信號(hào)碼元采樣M次，即Si=[si1,…,siM]T為M維向量。當(dāng)采用二進(jìn)制相移鍵控(BinaryPhaseShiftKeying,BPSK)調(diào)制時(shí)，sij∈{1,-1}。根據(jù)上述討論，接收模型可以表示為：

Xi=aSi+Ni(i=1,…,N)

(7)其中，a為信道衰減系數(shù)，M維向量Xi=[xi1,…,xiM]T為接收到的信號(hào)碼元，Ni=[ni1,…,niM]T為噪聲，根據(jù)假設(shè)③，在每次處理窗口期間a是固定的，噪聲的采樣值服從1.2節(jié)混合模型的分布，并滿足獨(dú)立同分布，則結(jié)合式(6)、式(7)，xij滿足：

(8)

2.2 貝葉斯層次模型

貝葉斯推斷通過(guò)先驗(yàn)分布來(lái)推斷后驗(yàn)分布：

(9)

根據(jù)貝葉斯層次理論，可以得出參數(shù)集變量{a,δ2,λ,σ2,ω,Z,S}的后驗(yàn)分布為：

(10)

對(duì)上述參數(shù)選取先驗(yàn)分布時(shí)，若參數(shù)存在共軛先驗(yàn)分布，則優(yōu)先選取共軛先驗(yàn)分布：

參數(shù)a的共軛先驗(yàn)分布為高斯分布：a～N(κ,ε2)。參數(shù)1/δ2的共軛先驗(yàn)分布為伽馬分布：1/δ2～Γ(α,β)。參數(shù)λ的先驗(yàn)分布為伽馬分布：λ～Γ(ξ,τ)。參數(shù)1/σ2的共軛先驗(yàn)分布為伽馬分布：1/σ2～Γ(g,h)。參數(shù)ω的共軛先驗(yàn)為對(duì)稱狄利克雷分布：ω～D(η,η)。對(duì)于信號(hào)碼元Si，通常認(rèn)為發(fā)射端各個(gè)信號(hào)碼元的發(fā)射概率相同，所以取先驗(yàn)分布：

(11)

為了直觀表述參數(shù)之間的關(guān)系，畫(huà)出模型參數(shù)的直接非循環(huán)圖，如圖1所示。

圖1 噪聲模型參數(shù)的直接非循環(huán)圖Fig.1 Direct acyclic graph specific to the parameters of the noise model

圖1中圓圈代表參量，是需要估計(jì)的值，方框代表定值。設(shè)θ={a,δ2,λ,σ2}，θ的超參數(shù)為φ={κ,ε2,α,β,ξ,τ,g,h}，則式(10)可以擴(kuò)展為：

(12)

3 盲檢測(cè)算法設(shè)計(jì)

盲檢測(cè)算法根據(jù)混合模型的性質(zhì)，采用MCMC算法實(shí)現(xiàn)信號(hào)碼元的盲檢測(cè)。在MCMC算法中，常采用Gibbs抽樣和Metropolis-Hasting(M-H)抽樣算法對(duì)更新參數(shù)進(jìn)行抽樣。本文中，對(duì)于存在共軛先驗(yàn)分布的參數(shù)采用Gibbs抽樣算法，對(duì)不存在共軛先驗(yàn)分布的參數(shù)采用M-H抽樣算法。算法在t步的流程如下：

步驟1：通過(guò)Gibbs抽樣更新參數(shù)ω。

ω的后驗(yàn)分布仍服從狄利克雷分布：

(13)

步驟2：通過(guò)Gibbs抽樣更新參數(shù)a。

衰減系數(shù)a的后驗(yàn)分布仍服從高斯分布：

(14)

其中

(15)

(16)

通過(guò)Gibbs抽樣，生成新的高斯分布隨機(jī)數(shù)對(duì)a進(jìn)行更新。

步驟3：通過(guò)Gibbs抽樣更新參數(shù)δ2。

1/δ2的后驗(yàn)分布為：

(17)

通過(guò)Gibbs抽樣，生成新的伽馬分布隨機(jī)數(shù)對(duì)δ2進(jìn)行更新。

步驟4：通過(guò)Gibbs抽樣更新參數(shù)σ2。

1/σ2的后驗(yàn)分布為：

(18)

通過(guò)Gibbs抽樣，生成新的伽馬分布隨機(jī)數(shù)對(duì)σ2進(jìn)行更新。

步驟5：通過(guò)M-H算法更新參數(shù)λ。

λ的后驗(yàn)概率為：

(19)

因?yàn)棣说暮篁?yàn)概率復(fù)雜，不是常見(jiàn)的分布，所以通過(guò)M-H抽樣算法更新λ值，算法步驟如下：

(b)從(a)中抽樣λ′作為備選值；

(c)從均勻分布中生成隨機(jī)數(shù)u～U(0,1)；

(d)計(jì)算新值的接受概率R=min(1,A)，其中：

(20)

(e)若u≤R，則接受新值λ(t+1)=λ′，否則不接受新值，則λ(t+1)=λ(t)。

步驟6：更新標(biāo)簽變量Z。

zij的后驗(yàn)概率為：

(21)

(22)

步驟7：通過(guò)M-H算法更新變量{γij}。

由于觀測(cè)數(shù)據(jù)只有一部分屬于高斯尺度混合分布，所以僅更新對(duì)應(yīng)這部分的γij，γij的后驗(yàn)概率為：

(23)

(a)γij選用建議分布為：γij～I(xiàn)G(5,λ)；

(c)從均勻分布中生成隨機(jī)數(shù)u～U(0,1)；

(d)計(jì)算新值的接受概率：

(24)

步驟8：通過(guò)Gibbs采樣更新信號(hào)碼元Si。

信號(hào)碼元Si的后驗(yàn)概率為：

(25)

(26)

算法在t+1步時(shí)，采用t步更新過(guò)的模型參數(shù)，重復(fù)上述8個(gè)步驟進(jìn)行新一次的更新，直至各參數(shù)收斂。參數(shù)收斂后的每次采樣值就可以看作目標(biāo)分布的采樣值，這部分采樣值的平均值就可以作為參數(shù)的估計(jì)值。

4 試驗(yàn)及驗(yàn)證

4.1 模型驗(yàn)證

圖2為實(shí)測(cè)的大氣噪聲數(shù)據(jù)，用本文的算法參數(shù)估計(jì)部分對(duì)噪聲進(jìn)行參數(shù)估計(jì)，設(shè)置μ=10，估計(jì)出的參數(shù)值為δ2=26.71，λ=1.60，σ2=6.88，ω={0.04,0.96}。圖3為用實(shí)測(cè)的大氣噪聲數(shù)據(jù)和估計(jì)參數(shù)畫(huà)出的幅度概率分布(AmplitudeProbabilityDistribution，APD)曲線，從圖中可以看出，估計(jì)參數(shù)畫(huà)出的APD曲線與實(shí)測(cè)的噪聲數(shù)據(jù)吻合的很好，平均相對(duì)誤差小于0.2%，說(shuō)明了本文的模型對(duì)接收機(jī)中實(shí)際噪聲的適用性。

圖2 實(shí)測(cè)大氣噪聲數(shù)據(jù)Fig.2 The actual noise data of the receiver

圖3 實(shí)測(cè)噪聲數(shù)據(jù)和估計(jì)參數(shù)下噪聲的APD曲線Fig.3 Amplitude probability distribution graph of the actual noise and noise with estimation parameters

4.2 參數(shù)估計(jì)

將混合模型的參數(shù)設(shè)置為δ2=2，λ=2，σ2=0.5，ω={0.2,0.8}，a=0.6,仿真數(shù)據(jù)數(shù)目設(shè)為N=2000，每個(gè)信號(hào)碼元采樣次數(shù)M=10。超參數(shù)簡(jiǎn)單地設(shè)置為κ=1，ε2=1，α=0.5，β=1，ξ=2，τ=1，g=0.5，h=1，b=0.5，η=0.5。設(shè)置算法的迭代次數(shù)為800，前200次為廢棄數(shù)據(jù)長(zhǎng)度。仿真的結(jié)果如圖4～8所示，因?yàn)棣?=1-ω2，所以只給出ω2的迭代收斂圖。

圖4 δ2的迭代收斂圖Fig.4 Time-average convergence of δ2

圖5 λ的迭代收斂圖Fig.5 Time-average convergence of λ

圖6 σ2的迭代收斂圖Fig.6 Time-average convergence of σ2

圖7 ω2的迭代收斂圖Fig.7 Time-average convergence of ω2

圖8 a的迭代收斂圖Fig.8 Time-average convergence of a

從各參數(shù)的迭代收斂情況可以看出，在進(jìn)行30次迭代后，各參數(shù)就收斂到實(shí)際值附近。對(duì)后600次的迭代值進(jìn)行平均，得到各參數(shù)的估計(jì)結(jié)果為δ2=2.18，λ=2.02，σ2=0.50，ω={0.19,0.81}，a=0.60。由于更新參數(shù)δ2和λ的過(guò)程中，采用了M-H算法，所以其迭代數(shù)據(jù)方差大，但均值仍然與實(shí)際值相符合。

4.3 誤碼率測(cè)試

為了測(cè)試盲檢測(cè)算法的性能，對(duì)最優(yōu)檢測(cè)和盲檢測(cè)算法的誤碼率進(jìn)行比較。最優(yōu)檢測(cè)為已預(yù)知模型的準(zhǔn)確參數(shù)，然后對(duì)信號(hào)碼元進(jìn)行檢測(cè)。基于逆高斯分布的高斯尺度混合模型的二階矩為μδ2，高斯分布的方差為σ2，所以總的噪聲功率為ω1μδ2+ω2σ2，由此噪聲功率來(lái)定義信噪比。

圖9為不同ω值下，盲檢測(cè)算法和最優(yōu)檢測(cè)的誤碼率性能比較。從圖中可以看出：①在ω值固定的條件下，盲檢測(cè)算法的性能在較高信噪比下都能逼近最優(yōu)檢測(cè)的性能，隨著信噪比的降低，較最優(yōu)檢測(cè)性能略有下降。②對(duì)于相同的信噪比，ω2的值越小，盲檢測(cè)算法的性能也隨著降低，這是因?yàn)?，?的值越小，即非高斯噪聲部分的能量越大，使得局部時(shí)段信噪比降低很多，從而導(dǎo)致誤碼率提高。③在信噪比很低的情況下，由于整個(gè)時(shí)段的信噪比都很低，所以ω2值不同對(duì)算法性能的影響不是很明顯。

圖9 不同ω值下盲檢測(cè)算法的誤碼率性能Fig.9 Performance of the blind detection algorithm with different parameter ω

圖10為不同過(guò)采樣率M下，盲檢測(cè)算法和最優(yōu)檢測(cè)的誤碼率性能比較。如圖10所示：無(wú)論過(guò)采樣率M取值為多少，盲檢測(cè)算法的誤碼率性能在較高信噪比下都逼近最優(yōu)檢測(cè)的誤碼率性能，隨著信噪比的降低，較最優(yōu)檢測(cè)性能略有下降；隨著M值變大，盲檢測(cè)算法的誤碼率變低，并且低信噪比情況下比高信噪比情況下，M值對(duì)誤碼率的影響較大。

圖10 不同過(guò)采樣率下盲檢測(cè)算法的誤碼率性能Fig.10 Performance of the blind detection algorithm with different oversampling rates

5 結(jié)論

根據(jù)實(shí)際超低頻接收機(jī)中噪聲的特點(diǎn)，對(duì)噪聲采用GSM分布和高斯分布的混合模型建模，在貝葉斯層次模型下，采用MCMC算法，通過(guò)Gibbs抽樣和M-H抽樣，同步估計(jì)信道衰落系數(shù)和噪聲模型參數(shù)，并對(duì)信號(hào)碼元進(jìn)行檢測(cè)。通過(guò)對(duì)實(shí)測(cè)噪聲數(shù)據(jù)分析，說(shuō)明該模型對(duì)接收機(jī)中實(shí)際接收的噪聲有很好的適用性。對(duì)算法性能進(jìn)行仿真分析，結(jié)果表明，MCMC算法迭代效率和精度高。在不同噪聲模型參數(shù)ω和過(guò)采樣率M下，盲檢測(cè)算法的性能在高信噪比下都逼近最優(yōu)檢測(cè)的性能，對(duì)超低頻非高斯噪聲下信號(hào)接收有實(shí)際的意義。本文算法有以下創(chuàng)新點(diǎn)：①對(duì)噪聲采用GSM分布和高斯分布的混合模型建模，符合實(shí)際接收機(jī)中的噪聲特點(diǎn)，對(duì)實(shí)際裝備應(yīng)用有現(xiàn)實(shí)意義和優(yōu)勢(shì)；②充分利用GSM分布的性質(zhì)，將GSM分布等價(jià)為條件高斯分布設(shè)計(jì)MCMC算法，使算法的迭代效率快、精度高。

References)

[1]ConteE，DiBisceglieM，LopsM．Optimumdetectionoffadingsignalsinimpulsivenoise[J]．IEEETransactionsonCommunications，1995，43(2/3/4)：869-876．

[2]BuzziS，ConteE，LopsM．Optimumdetectionoverrayleigh-fading，dispersivechannels，withnon-Gaussiannoise[J]．

IEEETransactionsonCommunications，1997，45(9)：1061-1069．

[3]BuzziS，ConteE，LopsM．Signaldetectionoverrayleigh-fadingchannelswithnon-Gaussiannoise[J]．IEEProceedingsCommunications，1997，144(6)：381-386．

[4] 應(yīng)文威, 蔣宇中, 劉月亮. 大氣低頻噪聲混合模型的MCMC參數(shù)估計(jì)[J]. 系統(tǒng)工程與電子技術(shù), 2012, 34(6): 1241-1245.

YINGWenwei,JIANGYuzhong,LIUYueliang.ParametersestimationformixturemodelofatmosphericnoisethroughMCMCmethod[J].SystemsEngineeringElectronics, 2012, 34(6): 1241-1245. (inChinese)

[5]YingWW,JiangYZ,LiuYL,etal.Anoveladaptivereceiverforflatfadingchannelswithimpulsivenoise[C]//ProceedingsofIEEE12thInternationalConferenceonComputerInformationTechnology, 2012: 631-634.

[6]OllilaE,TylerDE,KoivunenV,etal.Compound-Gaussiancluttermodelingwithaninversegaussiantexturedistribution[J].IEEESignalProcessingLetters, 2012, 19(12): 876-879.

[7]GrecoM,GiniF,DianiM.RobustCFARdetectionofrandomsignalsincompound-Gaussianclutterplusthermalnoise[J].IEEProceedingsRadar,SonarNavigation, 2001, 148(4): 227-232.

[8]WangJ,DogandzicA,NehoraiA.Maximumlikelihoodestimationofcompound-Gaussiancluttertargetparameters[J].IEEETransactionsonSignalProcessing, 2006, 54(10): 3884-3898.

[9] 趙宜楠, 李風(fēng)從, 尹彬. 嚴(yán)重拖尾復(fù)合高斯雜波中目標(biāo)的自適應(yīng)極化檢測(cè)[J]. 電子與信息學(xué)報(bào), 2013, 35(2): 376-380.

ZHAOYinan,LIFengcong,YINBin.Adaptivepolarimetricdetectionoftargetsinheavy-tailedcompound-Gaussianclutter[J].JournalofElectronics&InformationTechnology, 2013, 35(2): 376-380. (inChinese)

[10]BalleriA,NehoraiA,WangJ.Maximumlikelihoodestimationforcompound-Gaussianclutterwithinversegammatexture[J].IEEETransactionsonAerospaceElectronicSystems, 2007, 43(2): 775-779.

[11]GengP,HuangZT,WangFH,etal.SinglechannelblindsignalseparationwithBayesian-MCMC[C]//ProceedingsofInternationalConferenceonWirelessCommunications&SignalProcessing,IEEE, 2009: 1-4.

[12]LeeA,YauC,GilesMB,etal.OntheutilityofgraphicscardstoperformmassivelyparallelsimulationofadvancedMonteCarlomethods[J]JournalofComputationalGraphicalStatistics, 2010, 19(4): 769-789.

Symbol detection algorithm in non-Gaussian noise using Markov chain Monte Carlo method

FENG Shimin1， ZHOU Suihua1， YING Wenwei2

Consideringthatthereceiverwasnotonlyaffectedbythenon-GaussiannoisebutalsoaffectedbyitsinternalandexternalenvironmentofGaussiannoise,amixedmodelcomposedbynon-GaussiandistributionplusGaussiandistributionwasproposed.AblinddetectionalgorithmbasedonMarkovChainMonteCarlomethodwasdesignedaccordingtothepropertiesofthemixedmodel.Theblinddetectionalgorithmcouldestimatethechannelfadingcoefficient,parametersofnoisemodelandcoulddetectsignalelement.DetectsignalsbasedonBayesianhierarchicalmodelwasusingGibbssampleandM-Hsampleforparameterupdating.Thealgorithmhasahighiterativeefficiencyandprecision.Resultsshowthattheproposedblinddetectionalgorithmperformsaswellastheoptimaldetectionalgorithmandhasimportantrealisticsignificanceinsuperlow-freguencysignalreception.

non-Gaussiannoise;blinddetection;Gaussianscalemixture;mixedmodel

2014-08-29

國(guó)家自然科學(xué)基金資助項(xiàng)目(51109215)

馮士民(1986—)，男，河北石家莊人，博士研究生，E-mail：fengshimin_86@126.com；周穗華(通信作者)，男，教授，博士，博士生導(dǎo)師，E-mail：zhousuihua@hotmail.com

10.11887/j.cn.201504019

http://journal.nudt.edu.cn

TN

A