一種基于混合MCMC的時域波達(dá)方向估計算法*

李國彬，宋曉鷗，鈕嘉銘

(武警工程大學(xué) a.研究生大隊；b.信息工程學(xué)院，西安 710086)

0 引 言

對于寬帶陣列信號波達(dá)方向的估計，一般都是通過構(gòu)建寬帶信號的頻域模型，將問題轉(zhuǎn)換到頻域進(jìn)行處理，例如ISM類算法[1-2]和CSM類算法[3],但這兩類算法在運(yùn)行的過程中都會引入轉(zhuǎn)換誤差，運(yùn)算復(fù)雜度十分高，同時在處理一些寬帶和窄帶同時存在的混合信號上的表現(xiàn)也欠佳。如何估計寬帶和窄帶并存的復(fù)雜信號，降低模型的估計誤差，是目前波達(dá)方向估計的主要研究方向。2005年，William等人[4]提出了一種寬帶、窄帶同時適用的的時域信號模型，為波達(dá)方向估計提供了新的思路，但是該方法是基于低通信號重構(gòu)理論下建立的模型，對于信號的采樣頻率要求較高，工程上難以實現(xiàn)。2011年，趙擁軍等人[5]在此模型的基礎(chǔ)上進(jìn)行了改進(jìn)，提出了一種以帶通信號重構(gòu)理論為基礎(chǔ)的陣列信號處理時域模型，有效地降低了采樣頻率，同時對模型的適用性也進(jìn)行了改進(jìn)。這兩種波達(dá)方向估計的時域模型都是通過建立關(guān)于未知參數(shù)的貝葉斯后驗分布函數(shù)來進(jìn)行參數(shù)的估計，因此貝葉斯后驗分布函數(shù)的求解精度直接影響著參數(shù)估計的性能，而時域模型的后驗分布函數(shù)是復(fù)雜的，在求解的過程中必然會遇到多維積分的計算，這極大地增大了求解的難度。采用最大似然法進(jìn)行多維搜索雖然可以近似地求解出貝葉斯后驗分布函數(shù)的極大值，但這顯然是復(fù)雜的，不僅運(yùn)算量大而且求解的精度也難以滿足實際的需求。

近年來，由于馬爾科夫鏈蒙特卡洛(Markov Chain Monte Carlo,MCMC)方法在統(tǒng)計物理學(xué)上的出色表現(xiàn)，引起了越來越多通信學(xué)者的注意，因此也被應(yīng)用到信號的參數(shù)估計上來[6-8]。MCMC方法可以產(chǎn)生一個平穩(wěn)分布為任意函數(shù)的馬爾科夫鏈，基于這一點可以將MCMC方法運(yùn)用到貝葉斯后驗分布函數(shù)的求解上來，提高貝葉斯函數(shù)求解的效率和精度。但是馬爾科夫鏈的構(gòu)造過程并不是簡單的，馬爾科夫鏈的初始狀態(tài)以及提議函數(shù)的選擇都會影響馬爾科夫鏈的收斂速度速度和收斂效果。考慮到這一點，本文主要對MCMC方法中馬爾科夫鏈的構(gòu)造過程進(jìn)行了改進(jìn)：采用多鏈并行的方式進(jìn)行馬爾科夫鏈的構(gòu)造，將原本的長鏈劃分為多個短鏈，并且就此模型對狀態(tài)空間進(jìn)行了優(yōu)化；在生成馬爾可夫鏈的過程中，基于交叉的原則動態(tài)地選擇提議函數(shù)，且每條短鏈采用不同的交叉原則；由提議函數(shù)生成候選狀態(tài)的過程采用單元素和多元素混合移動的方式。

1 時域信號模型及貝葉斯后驗分布函數(shù)

為了便于分析，選擇陣列結(jié)構(gòu)較為簡單的M元均勻直線陣作為信號的接收陣列。假設(shè)有K個信號源從遠(yuǎn)場入射，第k個信號源的頻率范圍為fk∈[fk1,fk2],其中，k=1,2,…,K,fk2=fk1+Bk,fk1和fk2分別代表上限頻率和下限頻率,Bk為信號帶寬,信號源的入射角度θ和陣元間的延遲τ可以表示為

θ=[θ1,θ2,…,θK]T,

(1)

τ=[τ1,τ2,…,τK]T。

(2)

式中：

(3)

式中：d表示均勻直線陣的陣元間距，c表示光傳播速度。

根據(jù)帶通信號的重構(gòu)定理，第k個信號源發(fā)射的信號sk(t)可以表示為

(4)

(5)

式中：ψ(t)為重構(gòu)函數(shù),[-L,L]為實際進(jìn)行重構(gòu)的區(qū)間。為了保證重構(gòu)信號不混疊，采樣周期Ts滿足[9]

(6)

由式(4)便可表示出第m個陣元接收到第k個信號為

(7)

定義

hl(t)=ψ(lTs-t)ω(lTs-t)，

(8)

ω(l-t)為窗函數(shù)(這里選擇海明窗)，可以降低由于取近似值帶來的誤差。

對式(7)進(jìn)行時域離散采樣，并用n-l代替l可得

(9)

式中：sk(n-mτk)?sk(nTs-mτk)，hl(mτk)?ψ(lTs-mτk)ω(lTs-mτk)。則陣元m接收到的信號為

m=0，1，…，M-1。

(10)

式中：σω表示噪聲的方差，ωnoise(n)表示均值為0、方差為1的高斯白噪聲。

天線陣的總輸出y(n)=[y0(n),y1(n),…,yM-1(n)]T可表示為

(11)

式中：

(12)

a(n)=[s1(n),s2(n),…,sk(n)]T。

(13)

陣列接收信號的時域模型為

H0(τ)a(n)+σωωnoise(n)。

(14)

利用貝葉斯準(zhǔn)則來估計波達(dá)方向的核心就是建立關(guān)于未知參數(shù)的后驗分布函數(shù)。式(14)共有τ、K、σω、a(n)四個位置參數(shù)，分別代表信號延時、信源個數(shù)、噪聲參數(shù)、信號矢量，本文只考慮在信號源個數(shù)已知情況下的信號時延的估計，基于陣列接收信號的時域模型，結(jié)合未知參數(shù)的先驗分布，通過積分的方式消除無關(guān)參數(shù)，得到關(guān)于信號時延的貝葉斯后驗概率函數(shù)[4]：

(15)

式中：γ0和ν0為超參數(shù)，一般設(shè)置為0；N為快拍數(shù)；δ為與信噪比有關(guān)的參數(shù)；

(16)

(17)

π(τ|Z)便是關(guān)于信號時延的貝葉斯后驗分布函數(shù)，如果能求解出其最大值，那么便可以得到基于樣本統(tǒng)計量的信號延時的最佳估計，但這顯然不是容易的：關(guān)于信號延時的貝葉斯后驗概率密度函數(shù)是一個多峰值函數(shù)，利用譜峰搜索技術(shù)求解最大值勢必會帶來龐大的計算量，求解的精度也會受限于搜索步長。

2 基于混合MCMC的時域波達(dá)方向估計算法

MCMC方法最關(guān)鍵的就是轉(zhuǎn)移核的構(gòu)造，選取的轉(zhuǎn)移核應(yīng)當(dāng)使得構(gòu)造的馬爾科夫鏈滿足細(xì)節(jié)平穩(wěn)。Metropolis-Hastings(M-H)抽樣法是MCMC中常用的一種方法，通過選取合適的提議函數(shù)來抽取樣本，并以一定的的概率選擇是否接受候選狀態(tài)，直到樣本最后收斂。

圖1為M-H抽樣法產(chǎn)生的馬爾科夫鏈的流程圖，初始狀態(tài)和提議函數(shù)相當(dāng)于M-H抽樣法的輸入，初始狀態(tài)的不同會導(dǎo)致馬爾科夫鏈達(dá)到平穩(wěn)所需要的時間不同，提議函數(shù)的選擇則直接關(guān)系到抽樣結(jié)果，選擇遍歷性不同的的提議函數(shù)會產(chǎn)生不同的收斂速率和結(jié)果。因此，我們考慮從初始狀態(tài)和提議函數(shù)入手來提高馬爾科夫鏈的收斂效果和速率。

圖1 M-H方法的流程圖

2.1 并行馬爾科夫鏈

以信源個數(shù)等于3的情況為例，結(jié)合式(3)我們可得出需要構(gòu)造的馬爾科夫鏈的狀態(tài)空間為一個“正方體”，在整個“正方體”中進(jìn)行馬爾科夫鏈樣本的選擇不僅需要的時間長，而且很難確保在有限迭代次數(shù)下構(gòu)造的馬爾科夫鏈能夠收斂。初始狀態(tài)的選擇不當(dāng)也可能造成馬爾科夫鏈未能收斂，雖然可以通過預(yù)估計的方式對初始狀態(tài)進(jìn)行修正，但這顯然是復(fù)雜的，對估計的實時性也會造成影響?？紤]到以上的因素，我們將馬爾科夫鏈的構(gòu)造過程拆分成多個并行的短鏈來進(jìn)行，也就是說將狀態(tài)空間劃分為多個子空間如圖2所示，在每個子空間分別進(jìn)行馬爾可夫鏈的構(gòu)造，構(gòu)造過程中的初始狀態(tài)是由每個子空間決定的，最后將每個子空間的結(jié)果進(jìn)行合并形成最終的馬爾科夫鏈。這種并行構(gòu)造方式極大地改善了初始狀態(tài)選擇不佳對馬爾科夫鏈?zhǔn)諗克俣群徒Y(jié)果的影響，在相同的迭代次數(shù)下馬爾科夫鏈也更容易達(dá)到收斂狀態(tài)。

圖2 并行構(gòu)造方式的狀態(tài)空間

基于此類模型，我們還發(fā)現(xiàn)貝葉斯后驗概率函數(shù)對于狀態(tài)空間的選擇是“無序”的。例如，實際時延值為[τ1,τ2,τ3]，而狀態(tài)空間中的其他“有序”狀態(tài)[τ1,τ3,τ2]、[τ2,τ1,τ3]、[τ2,τ3,τ1]等，與實際時延[τ1,τ2,τ3]有著相同的貝葉斯后驗概率函數(shù)值，這對馬爾科夫鏈的構(gòu)造是非常不利的，嚴(yán)重影響了算法的收斂速率。因此，我們采用圖3的方法對狀態(tài)空間進(jìn)行了優(yōu)化和改進(jìn)，可以看出改進(jìn)后的狀態(tài)空間相比于之前明顯減少，但是又包含了所有的“無序”狀態(tài)。事實上，經(jīng)過改進(jìn)后的狀態(tài)空間依然存在部分“有序”狀態(tài)，但去除這部分“有序”狀態(tài)會使算法的復(fù)雜度驟升，因此我們選擇保留這部分“有序”狀態(tài)，并且少量的“有序”狀態(tài)對于算法的收斂速度和收斂結(jié)果也有一定的提升。

圖3 改進(jìn)的狀態(tài)空間圖

2.2 混合MCMC方法

選擇遍歷性強(qiáng)弱不同的提議函數(shù)對馬爾科夫鏈的收斂速度和結(jié)果影響十分大，遍歷性弱的提議函數(shù)會導(dǎo)致算法收斂速度慢，遍歷性強(qiáng)的提議函數(shù)又可能會出現(xiàn)近似收斂的情況?？紤]到這點，這里采用一種混合的提議函數(shù)選取方式：在生成馬爾科夫鏈的過程中，隨機(jī)交叉地使用遍歷性弱的隨機(jī)游走采樣(下一次的采樣結(jié)果是當(dāng)前狀態(tài)的隨機(jī)擾動)和遍歷性強(qiáng)的獨立游走采樣(下一次的采樣結(jié)果是獨立的，與當(dāng)前狀態(tài)無關(guān))。這種混合的采樣方式相比于隨機(jī)游走采樣收斂速度更快，相比于獨立游走采樣收斂結(jié)果更加穩(wěn)定。從圖4可以看出，由提議函數(shù)生成候選狀態(tài)的過程中，提議函數(shù)的選擇是隨機(jī)多變的，包含了遍歷性強(qiáng)弱不同的提議函數(shù)。

圖4 混合MCMC的算法流程圖

同時此模型下的候選狀態(tài)是一個多維向量，每個維度的值對該狀態(tài)都有貢獻(xiàn)，如果我們同時生成多維度的候選狀態(tài)，就很容易把一些單維度性能較好的狀態(tài)給拒絕，而單元素移動每次只改變一個維度的值，這樣就會將收斂性較好的維度傳遞下去，從而加快收斂速度和改善收斂結(jié)果[10]。改進(jìn)后的算法采用多元素和單元素混合移動的辦法來生成候選狀態(tài)，如圖4，在馬爾科夫鏈未達(dá)到收斂狀態(tài)前采用多元素移動的方式來生成候選狀態(tài)，加速狀態(tài)的游離，在馬爾科夫鏈接近收斂狀態(tài)后(一般根據(jù)經(jīng)驗來決定)采用單元素。

3 仿真與分析

實驗采用陣元數(shù)目M=12的均勻直線陣，陣元間距為波長的一半，時域模型中的插值個數(shù)2L-1=11，信號源的個數(shù)K=3(其中兩個為窄帶調(diào)幅信號，一個為寬帶線性調(diào)頻信號,信號的強(qiáng)度比為2∶3∶4)，信號源的方向分別為10°、15°、20°。采樣的快拍數(shù)為128，仿真環(huán)境中的信噪比為10 dB，馬爾科夫鏈的迭代總次數(shù)為2 000次。

圖5給出了采用單鏈構(gòu)造方式的估計結(jié)果，可以看出馬爾科夫鏈需要大約1 800次的迭代樣本才能夠收斂。圖6給出了采用多鏈并行方式的估計結(jié)果，可見在實際波達(dá)方向?qū)?yīng)的子空間2內(nèi)，只需要450次左右迭代樣本便可以收斂。事實上采用并行的方式可以在相同的時間完成多個(短)鏈的構(gòu)造，因此算法的實時性并不會降低，算法的收斂結(jié)果反而會更加精確。并且我們可以看出，雖然實際的波達(dá)方向并不其余在子空間內(nèi)，但是其他子空間構(gòu)造的馬爾科夫鏈仍然有部分正確結(jié)果，也就是說單維度的收斂性也能夠傳遞下來，這就為單元素移動提供了一定的理論依據(jù)。

圖5 MCMC方法的估計結(jié)果

圖6 多鏈并行方式下的估計結(jié)果

采用混合MCMC的方式構(gòu)造馬爾科夫鏈，在1 800次迭代后樣本便可達(dá)到收斂。作為對比，圖7給出了采用遍歷性強(qiáng)的獨立游走采樣的估計結(jié)果，可以看出經(jīng)過2 000次迭代后，樣本達(dá)到了一種“偽收斂”的狀態(tài)，這對于波達(dá)方向的估計來講是非常不利的。

圖7 獨立游走采樣的估計結(jié)果

圖8給出了在迭代次數(shù)末端分別使用單元素移動和多元素移動的DOA估計結(jié)果，可以看出，在馬爾科夫鏈達(dá)到收斂狀態(tài)后單元素移動的估計結(jié)果在真實值附近波動的幅度更小，這也是為什么在馬爾科夫鏈的末端使用單元素移動的原因。

圖8 使用多元素移動和單元素移動的DOA估計結(jié)果

表1給出了1 000次蒙特卡洛實驗中單元素移動和多元素移動的估計結(jié)果和均方誤差(Mean Squared Error，MSE)[11],可以看出單元素移動的方式極大地改善了MCMC方法的估計精確度。

表1 估計結(jié)果和誤差

(18)

圖9給出了不同信噪下和快拍數(shù)下本文算法與文獻(xiàn)[5]算法和MUSIC算法DOA估計結(jié)果的均方誤差，可以看出當(dāng)快拍數(shù)固定為128時，MUSIC算法在信噪比較低的情況下表現(xiàn)出極差的性能，而文獻(xiàn)[5]算法和改進(jìn)的并行混合算法在不同信噪比下都表現(xiàn)出較好性能，估計的均方誤差遠(yuǎn)小于MUSIC算法，同時改進(jìn)后的并行算法相比于文獻(xiàn)[5]算法的估計性能又有了一定的提升；當(dāng)信噪比固定為10 dB時，MUSIC算法的估計性能隨快拍數(shù)的增加并沒有發(fā)生明顯的變化，而文獻(xiàn)[5]算法和改進(jìn)的并行混合算法都發(fā)生了明顯的變化，這與其算法的本質(zhì)有關(guān)(文獻(xiàn)[5]算法和本文的算法本質(zhì)上都是利用了貝葉斯原理求解，MUSIC算法則是利用信號與噪聲的相關(guān)特性)，并且改進(jìn)的并行混合算法的估計性能要略好于文獻(xiàn)[5]算法。

4 結(jié)束語

通過構(gòu)造馬爾科夫鏈的方法可以實現(xiàn)高分辨率的波達(dá)方向估計，但是馬爾科夫鏈對初始狀態(tài)、提議函數(shù)等有很高的要求，因此本文采用了一種并行混合的方式進(jìn)行構(gòu)造。通過實驗仿真可知，本文的模型可同時適用于寬帶和窄帶信號的波達(dá)方向估計；在相同的迭代次數(shù)下，本文的并行混合構(gòu)造方式相比于文獻(xiàn)[5]算法表現(xiàn)出更好的收斂穩(wěn)定性，估計的精度也要高于傳統(tǒng)的MUSIC算法和文獻(xiàn)[5]的算法；同時，本文采用的并行運(yùn)算方式很大程度上減少了算法的運(yùn)算時間，提高了估計的實時性。此算法可用于復(fù)雜信號的高分辨率波達(dá)方向估計，估計的精度和實時性都能滿足實際需求。

下一步工作將研究基于波達(dá)方向估計的自適應(yīng)抗干擾技術(shù)，將此算法的估計結(jié)果與自適應(yīng)抗干擾技術(shù)相結(jié)合，提高系統(tǒng)的抗干擾能力。