馮士民，周穗華，應文威

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

馮士民1，周穗華1，應文威2

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

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

非高斯噪聲；盲檢測；高斯尺度混合；混合模型

(1.DepartmentofWeaponryEngineering,NavalUniversityofEngineering,Wuhan430033,China;

2.ThePLAUnitof91635,Beijing102249,China)

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

1 噪聲模型

1.1 GSM模型

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

(1)

n～N(0,γδ2), γ～IG(μ,λ)

(2)

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

(3)

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

1.2 噪聲模型

為了對噪聲進行準確估計并進行信號碼元檢測，采用高斯尺度混合分布和高斯分布的混合模型作為噪聲模型：

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

(4)

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

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

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

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

(5)

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

(6)

2 信號碼元檢測模型

2.1 接收信號模型

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

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

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

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

(8)

2.2 貝葉斯層次模型

貝葉斯推斷通過先驗分布來推斷后驗分布：

(9)

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

(10)

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

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

(11)

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

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

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

(12)

3 盲檢測算法設計

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

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

ω的后驗分布仍服從狄利克雷分布：

(13)

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

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

(14)

其中

(15)

(16)

通過Gibbs抽樣，生成新的高斯分布隨機數(shù)對a進行更新。

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

1/δ2的后驗分布為：

(17)

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

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

1/σ2的后驗分布為：

(18)

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

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

λ的后驗概率為：

(19)

因為λ的后驗概率復雜，不是常見的分布，所以通過M-H抽樣算法更新λ值，算法步驟如下：

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

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

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

(20)

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

步驟6：更新標簽變量Z。

zij的后驗概率為：

(21)

(22)

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

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

(23)

(a)γij選用建議分布為：γij～IG(5,λ)；

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

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

(24)

步驟8：通過Gibbs采樣更新信號碼元Si。

信號碼元Si的后驗概率為：

(25)

(26)

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

4 試驗及驗證

4.1 模型驗證

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

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

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

4.2 參數(shù)估計

將混合模型的參數(shù)設置為δ2=2，λ=2，σ2=0.5，ω={0.2,0.8}，a=0.6,仿真數(shù)據(jù)數(shù)目設為N=2000，每個信號碼元采樣次數(shù)M=10。超參數(shù)簡單地設置為κ=1，ε2=1，α=0.5，β=1，ξ=2，τ=1，g=0.5，h=1，b=0.5，η=0.5。設置算法的迭代次數(shù)為800，前200次為廢棄數(shù)據(jù)長度。仿真的結(jié)果如圖4～8所示，因為ω1=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ù)的迭代收斂情況可以看出，在進行30次迭代后，各參數(shù)就收斂到實際值附近。對后600次的迭代值進行平均，得到各參數(shù)的估計結(jié)果為δ2=2.18，λ=2.02，σ2=0.50，ω={0.19,0.81}，a=0.60。由于更新參數(shù)δ2和λ的過程中，采用了M-H算法，所以其迭代數(shù)據(jù)方差大，但均值仍然與實際值相符合。

4.3 誤碼率測試

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

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

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

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

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

5 結(jié)論

根據(jù)實際超低頻接收機中噪聲的特點，對噪聲采用GSM分布和高斯分布的混合模型建模，在貝葉斯層次模型下，采用MCMC算法，通過Gibbs抽樣和M-H抽樣，同步估計信道衰落系數(shù)和噪聲模型參數(shù)，并對信號碼元進行檢測。通過對實測噪聲數(shù)據(jù)分析，說明該模型對接收機中實際接收的噪聲有很好的適用性。對算法性能進行仿真分析，結(jié)果表明，MCMC算法迭代效率和精度高。在不同噪聲模型參數(shù)ω和過采樣率M下，盲檢測算法的性能在高信噪比下都逼近最優(yōu)檢測的性能，對超低頻非高斯噪聲下信號接收有實際的意義。本文算法有以下創(chuàng)新點：①對噪聲采用GSM分布和高斯分布的混合模型建模，符合實際接收機中的噪聲特點，對實際裝備應用有現(xiàn)實意義和優(yōu)勢；②充分利用GSM分布的性質(zhì)，將GSM分布等價為條件高斯分布設計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] 應文威, 蔣宇中, 劉月亮. 大氣低頻噪聲混合模型的MCMC參數(shù)估計[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] 趙宜楠, 李風從, 尹彬. 嚴重拖尾復合高斯雜波中目標的自適應極化檢測[J]. 電子與信息學報, 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

國家自然科學基金資助項目(51109215)

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

10.11887/j.cn.201504019

http://journal.nudt.edu.cn

TN

A