基于迂回式匹配追蹤算法的DOA估計(jì)方法

黃嘉俊，張靖奇，何偉嘉，王 鵬

(中北大學(xué) 理學(xué)院，山西 太原 030051)

0 引 言

水聲信號的波達(dá)方向(Direction of arrival，DOA)估計(jì)一直是水聲信號處理的重要研究領(lǐng)域.早期的DOA估計(jì)算法主要圍繞波束形成技術(shù)開展，以Capon算法[1]為代表的波束形成算法，通過對傳感器陣列接收到的各路信號進(jìn)行加權(quán)處理，抑制非目標(biāo)方向的干擾信號，增強(qiáng)目標(biāo)方向的期望信號來獲得空間譜.然而，該技術(shù)無法在低信噪比和高混響環(huán)境中準(zhǔn)確估計(jì)波達(dá)方向.隨著子空間類算法的出現(xiàn)，波達(dá)方向估計(jì)的分辨率得到了極大的提高.一類是子空間分解類算法，包括多重信號分類(Multiple signal classification,MUSIC)算法[2]和旋轉(zhuǎn)不變子空間(Estimation of signal parameters via rotation invariant technique,ESPRIT)算法[3].這類算法利用信號子空間與噪聲子空間的正交特性來構(gòu)造空間譜，突破了DOA估計(jì)中分辨率的瑞利限，實(shí)現(xiàn)了超分辨率的DOA估計(jì).但是子空間分解類算法由于依賴陣列接收信號的統(tǒng)計(jì)特性，需要大量的快拍數(shù)據(jù)，難以實(shí)現(xiàn)小快拍數(shù)據(jù)下精確的空間譜估計(jì)，此外，這類方法無法有效估計(jì)相干信號的空間譜.另一類是子空間擬合類DOA估計(jì)算法，其中具有代表性的是最大似然(Maximum likelihood，ML)算法[4]，該方法通過求解非線性方向估計(jì)似然函數(shù)來獲得空間譜，由于需要進(jìn)行多維搜索來找尋最優(yōu)解，其運(yùn)算量較大.

近年來，Donoho D、Candes E和Tao T提出了一種新理論—壓縮感知(Compressive sensing，CS)[5].該方法可以在確保信息完整的情況下，以小于奈奎斯特采樣率的頻率對信號進(jìn)行采樣.其過程主要是：首先利用信號的稀疏性，以遠(yuǎn)小于奈奎斯特采樣率的隨機(jī)采樣率對稀疏信號進(jìn)行采樣，然后通過重構(gòu)算法較好地重構(gòu)原始信號.在實(shí)際空域中，實(shí)際目標(biāo)信號數(shù)相對于整個空間潛在的輻射信號數(shù)來說是極少的，這使得DOA估計(jì)中的空域信號本身即滿足稀疏性.因此，基于壓縮感知原理的DOA估計(jì)算法是一種利用信號稀疏性來實(shí)現(xiàn)高分辨率估計(jì)的方法[6].由于其在單快拍以及信源信號相干的條件下仍能獲得良好的DOA估計(jì)而受到廣泛的關(guān)注.

學(xué)者們在此基礎(chǔ)上提出了許多信號重構(gòu)的方法，主要包括凸優(yōu)化算法以及貪婪算法[7-8]等.基于凸優(yōu)化類方法的信號重構(gòu)算法，包括奇異值分解(Singular value decomposition，SVD)[9]、基追蹤降噪(Basis pursuit de-noising，BPDN)[10]和硬閾值迭代(Iterative hard thresholding，IHT)[11]等算法.該類算法是一種全局優(yōu)化算法，通過將非凸問題轉(zhuǎn)化為凸問題求解找到待重構(gòu)信號的近似值.文獻(xiàn)[12]給出了一種基于同倫思想的基追蹤問題解法，該方法通過迭代得出最有效的規(guī)整因子和重構(gòu)結(jié)果，適合數(shù)據(jù)量較少、無法確定誤差界限的場景，聯(lián)系寬帶信號的聯(lián)合稀疏性，該方法被擴(kuò)展到單個頻域快拍時的寬帶信號測向中.文獻(xiàn)[13]提出使用基追蹤降噪算法進(jìn)行DOA估計(jì)，基于壓縮感知的DOA估計(jì)模型的求解為最小l0范數(shù)問題.在一定條件下，l1范數(shù)與l0范數(shù)最小化問題是等價的[5]，該算法的思路是將l0范數(shù)最小化問題轉(zhuǎn)換成l1范數(shù)最優(yōu)化問題，以便使用最優(yōu)化方法求解.基追蹤降噪算法可以對相干信號進(jìn)行DOA估計(jì)，具有較高的穩(wěn)定性和重構(gòu)精度，但是其最終求解的無約束優(yōu)化問題涉及到正則化參數(shù)的選擇，正則化參數(shù)的大小影響著最終重構(gòu)效果對噪聲的適應(yīng)性，需要選擇合適的正則化參數(shù)來保證算法性能.迂回式匹配追蹤(Detouring matching pursuit，DMP)算法[14]是由裴廷睿等提出的一種貪婪重構(gòu)算法.該算法嚴(yán)格證明了子內(nèi)積逆和系數(shù)矩陣遞增遞減核心式，利用子內(nèi)積逆和系數(shù)矩陣減少計(jì)算殘差誤差變化量的計(jì)算量，降低了計(jì)算復(fù)雜度，具有重構(gòu)準(zhǔn)確率高、計(jì)算復(fù)雜度低、對傳感矩陣列的相關(guān)性要求低等特點(diǎn).

綜合上述分析，為了實(shí)現(xiàn)矢量水聽器陣列在小快拍及低信噪比下的波達(dá)方向估計(jì)，本文將迂回式匹配追蹤算法與波達(dá)方向估計(jì)相結(jié)合，提出一種基于迂回式匹配追蹤的波達(dá)方向估計(jì)算法.在信號重構(gòu)過程中，采用先逐個最優(yōu)縮減、后逐個最優(yōu)擴(kuò)增的方式更新支撐集，提高重構(gòu)準(zhǔn)確率，擴(kuò)大重構(gòu)稀疏信號的稀疏度范圍.在存在相干信號、相同的低信噪比、小快拍數(shù)條件下的仿真實(shí)驗(yàn)中，本方法可以得到有效的DOA估計(jì)，與傳統(tǒng)超分辨率算法相比更有優(yōu)勢.在真實(shí)的水聲實(shí)驗(yàn)環(huán)境中，應(yīng)用本文方法仍可以得到有效的DOA估計(jì).

1 矢量水聽器陣列信號模型

考慮到水聲環(huán)境的特殊性，假設(shè)聲信號為遠(yuǎn)場窄帶信號；水下聲速為1 500 m/s；水下傳播介質(zhì)為各向同性的非色散均勻線性介質(zhì)；仿真實(shí)驗(yàn)中使用的噪聲為高斯白噪聲；陣列的各傳感器間不存在互耦情形.

假設(shè)t時刻有K個空間信號s1(t),s2(t),…,sK(t)入射到由M個矢量水聽器按間距d組成的均勻線陣中，第k個信號的到達(dá)角為θk,k=1,2,…,K.以陣列中的第一個矢量水聽器為參考，t時刻第m個陣元的輸出為

m=1,2,…,M,

(1)

(2)

(3)

式中：k=1,2,…,K，將該陣列的輸出寫成矩陣形式，即

X(t)=A(θ)S(t)+N(t),

(4)

式中：X(t)=[x1(t),x2(t),…,xM(t)]T，X(t)∈CM×1為陣列輸出向量；A(θ)=[a(θ1)?u1,a(θ2)?u2,…,a(θK)?uK]為矢量水聽器陣列的信號方向矩陣，也叫流形矩陣；符號?為克羅內(nèi)克積；S(t)=[s1(t),s2(t),…,sK(t)]T，S(t)∈CK×1稱為信號源矢量；N(t)=[n1(t),n2(t),…,nM(t)]T，N(t)∈CM×1為陣列接收噪聲矢量.

2 基于迂回式匹配追蹤算法的DOA估計(jì)

2.1 DOA估計(jì)模型的稀疏化描述

G(θ)=[a(θ1),a(θ2),…,a(θN)]，

(5)

故DOA估計(jì)的數(shù)學(xué)模型可以轉(zhuǎn)化為

X(t)=G(θ)S(t)+N(t).

(6)

該式為典型的稀疏表示模型，因?yàn)镾(t)為空域內(nèi)的稀疏信號，不需要額外的稀疏基來對S(t)稀疏表示，所以擴(kuò)展后的超完備冗余字典G(θ)既是壓縮感知模型中的測量矩陣也是傳感矩陣.同時，測量矩陣是通過對空間進(jìn)行等角度網(wǎng)格劃分得到的，由這種稀疏表示方式得到的測量矩陣滿足等距約束性質(zhì)(Restricted isometry property,RIP)，因此可以通過求解式(6)來解決稀疏信號重構(gòu)問題，從而得到DOA估計(jì)值.

然而，式(6)是一個欠定方程，無法直接求解，信號的重構(gòu)為求解l0范數(shù)最小化問題，即通過

s.t.X(t)=G(θ)S(t)+N(t),

(7)

進(jìn)行求解.

然而，l0范數(shù)最小化為NP-hard問題，無法直接求解，具有代表性的解法主要有凸優(yōu)化算法和貪婪類重構(gòu)算法.凸優(yōu)化算法是一種全局優(yōu)化算法，如基追蹤(BP)算法.凸優(yōu)化算法通常將式(7) 轉(zhuǎn)化為帶有不等式約束的優(yōu)化問題，即最小l1范數(shù)優(yōu)化問題

(8)

式中：γ為誤差.一般地，將式(8)轉(zhuǎn)化成

(9)

式中：λ是一個與噪聲有關(guān)的正則化參數(shù)，即求解無約束優(yōu)化問題.

然而，最小l1范數(shù)優(yōu)化問題是一個廣義線性規(guī)劃問題，需要求解二階錐規(guī)劃，以犧牲計(jì)算復(fù)雜度來實(shí)現(xiàn)準(zhǔn)確重構(gòu)，計(jì)算復(fù)雜度較高，難以適用于大規(guī)模信號的快速DOA估計(jì).

2.2 DMP算法

輸入：M×N維傳感矩陣G(θ)，M維觀測信號X(t)，稀疏信號維度N，稀疏度K，閾值ε.

步驟1(初始化)：

步驟2(獲取初始值支撐集)：

步驟3：

令b=b+1.判定是否滿足b

步驟5(擴(kuò)增支撐集)：

νt=

Δ(Λ/{ω1,ω2,…,ωb})∪{νb,…,νt}).

選擇νt加入支撐集Λ，判定{ωb,…,ωt}∪{νb,…,νt}和{ωb,…,ωt}所包含的元素是否相同，若不相同，則t=t-1，計(jì)算子內(nèi)積逆矩陣P(Λ/{ω1,…,ωb})∪{νb,…,νt}和系數(shù)矩陣C(Λ/{ω1,…,ωb})∪{νb,…,νt}，重復(fù)步驟5，直至t=0；若相同，則跳轉(zhuǎn)至步驟3.

步驟6：

步驟7(更新支撐集)：

2.3 計(jì)算DOA估計(jì)值

3 算法復(fù)雜度分析

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

4.1 DMP算法、MUSIC算法在小快拍數(shù)、低信噪比條件下的估計(jì)性能

在該實(shí)驗(yàn)中，分別設(shè)置單信源與多信源的環(huán)境.其中陣元數(shù)目為10，陣元間距為0.5 m，DMP算法為單快拍，MUSIC算法快拍數(shù)為10，信噪比(Signal-to-noise ratio，SNR)均為-10 dB，單信源實(shí)驗(yàn)的入射角為-15°，信號頻率為 1 500 Hz.多信源實(shí)驗(yàn)的3個入射角分別為-15°，30°，50°，信號頻率分別為1 000 Hz，1 500 Hz，2 000 Hz.實(shí)驗(yàn)結(jié)果如圖1 所示.

由圖1 可以看出，在快拍數(shù)為10，SNR為-10 dB 的環(huán)境中，無論是單信源還是多信源實(shí)驗(yàn)，MUSIC算法和ESPRIT算法已經(jīng)無法估計(jì)出正確的DOA值，MUSIC算法還產(chǎn)生了大量的偽峰.在圖1(a)中標(biāo)記極值點(diǎn)的位置也可以看出，MUSIC算法在單信源數(shù)下搜尋到的譜峰為-5°，ESPRIT算法得出的估計(jì)角度為-6.262 25°，與實(shí)際的DOA值誤差較大.而本文所提出的DMP算法在單快拍條件下，估計(jì)信號入射角度為-14.5°，與實(shí)際信號入射角誤差為0.5°，依舊準(zhǔn)確估計(jì)出了入射信號方位角.在圖1(b)中，MUSIC算法在多信源數(shù)下搜尋到的譜峰為-17.5°，8.5°和46°，ESPRIT算法得出的估計(jì)角度為-11.861 4°，44.388 8°和65.119 8°.DMP算法在單快拍條件下，估計(jì)信號入射角度為-15.5°，30.5°和48.5°，與實(shí)際信號入射角的誤差介于0.5°～1.5°之間，表現(xiàn)出較好的估計(jì)效果.

(a) 單信源實(shí)驗(yàn)

(a) 1個信號源

為了說明DMP算法在單快拍、低信噪比環(huán)境中的優(yōu)勢，分別將MUSIC算法和ESPRIT算法的快拍數(shù)取為10，20，50和100，其他條件與上述實(shí)驗(yàn)一樣，表1 展示了不同快拍數(shù)下MUSIC算法和ESPRIT算法的單信源實(shí)驗(yàn)DOA估計(jì)結(jié)果.

表1 不同快拍數(shù)下MUSIC算法和ESPRIT算法的單信源實(shí)驗(yàn)DOA估計(jì)結(jié)果Tab.1 Single source experiment DOA estimation results of MUSIC and ESPRIT under different snapshots

由表1 可知，當(dāng)快拍數(shù)增加到100時，MUSIC和ESPRIT算法均可估計(jì)出正確的DOA值，但是在實(shí)際應(yīng)用中，要得到持續(xù)不斷的高質(zhì)量采樣往往是困難的，DMP算法可以在小快拍情形下得到正確的DOA估計(jì)值，有效減小了MEMS矢量水聽器陣列在數(shù)據(jù)采樣時的壓力，說明DMP算法具有實(shí)際應(yīng)用的價值.

4.2 不同信噪比下各算法估計(jì)誤差的比較

設(shè)置矢量均勻線陣陣元數(shù)為10，快拍數(shù)為10，信噪比以5 dB為步長從-10 dB增加到20 dB，每個信噪比下進(jìn)行100次蒙特卡洛實(shí)驗(yàn).圖2 展示了1個信號源、2個信號源和3個信號源條件下，DMP算法、MUSIC算法和ESPRIT算法各自的均方根誤差隨信噪比變化的曲線.

由圖2 可知，隨著信噪比增大，3種算法的均方根誤差都在逐漸減小，估計(jì)精度都在逐漸增大.在SNR處于-10 dB～-5 dB區(qū)間時，DMP算法的估計(jì)性能與MUSIC和ESPRIT算法接近，但考慮到DMP算法只取了單快拍下的數(shù)據(jù)進(jìn)行DOA估計(jì)，因此在信噪比較低時DMP算法只需要一次采樣便可達(dá)到多快拍下MUSIC和ESPRIT算法的估計(jì)性能.-5 dB以上時DMP算法的性能明顯優(yōu)于MUSIC和ESPRIT算法.從整體來看，DMP算法的均方根誤差都處于比較低的位置，估計(jì)精度較高.總之，DMP算法具有良好的抗噪能力，相同信噪比下只需單快拍數(shù)據(jù)即可取得相對較高的估計(jì)精度，在信噪比較高的環(huán)境下具有更顯著的優(yōu)勢.

4.3 相干信號源條件下的DOA估計(jì)性能對比

設(shè)置矢量均勻線陣陣元數(shù)為10，間距為0.5 m，3個聲源信號的入射角分別為-50°，10°，60°，其中10°與60°方向信號為相干信號，頻率分別為1 000 Hz，1 500 Hz和1 500 Hz.為了觀察相干信號對MUSIC算法的影響，將MUSIC算法的快拍數(shù)設(shè)置為100，DMP算法為單快拍，信噪比為0 dB，實(shí)驗(yàn)結(jié)果見圖3.

圖3 相干信號的DOA估計(jì)實(shí)驗(yàn)Fig.3 DOA estimation experiment of coherent signal

由圖3 可以看出，當(dāng)存在相干信源時，MUSIC算法出現(xiàn)了一些偽峰，并且其譜峰較為平滑，幾乎無法分辨信源角度.ESPRIT算法則無法得出正確的DOA估計(jì)值.DMP算法在單快拍、低信噪比情形下仍能得到有效的DOA估計(jì)值，這是由于DMP算法僅僅利用信號源的稀疏性來求解DOA估計(jì)值，不受入射信號頻率的影響，所以相干信號不會對DOA估計(jì)的結(jié)果產(chǎn)生較大影響.

5 結(jié) 論

針對單快拍、信噪比較低條件下的波達(dá)方向估計(jì)問題，本文提出一種基于迂回式匹配追蹤算法的DOA估計(jì)算法.該算法的優(yōu)勢和特點(diǎn)為：

1) 利用壓縮感知理論求解DOA問題，通過單快拍即可求解出聲源信號的入射角度，尤其是在低信噪比情形下，單快拍即可達(dá)到小快拍數(shù)下傳統(tǒng)的子空間分解類方法的估計(jì)精度.

2) 仿真實(shí)驗(yàn)表明，基于迂回式匹配追蹤算法的DOA估計(jì)算法求解的過程不涉及矩陣的特征分解，克服了相干信號對估計(jì)精度造成的影響，在處理小塊拍、低信噪比的相干信號時仍然有效.