李國(guó)彬，宋曉鷗，鈕嘉銘

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

0 引 言

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

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

1 時(shí)域信號(hào)模型及貝葉斯后驗(yàn)分布函數(shù)

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

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

(1)

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

(2)

式中：

(3)

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

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

(4)

(5)

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

(6)

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

(7)

定義

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

(8)

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

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

(9)

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

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)

陣列接收信號(hào)的時(shí)域模型為

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

(14)

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

(15)

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

(16)

(17)

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

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

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

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

圖1 M-H方法的流程圖

2.1 并行馬爾科夫鏈

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

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

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

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

2.2 混合MCMC方法

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

圖4 混合MCMC的算法流程圖

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

3 仿真與分析

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

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

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

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

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

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

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

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

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

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

(18)

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

4 結(jié)束語(yǔ)

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

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