賈思宇，路 茗，丁華澤，陳 明，趙魯陽

（1.中國科學(xué)院上海微系統(tǒng)與信息技術(shù)研究所無線傳感網(wǎng)與通信重點(diǎn)實(shí)驗(yàn)室，上海 200050；2.上?？萍即髮W(xué)信息科學(xué)與技術(shù)學(xué)院，上海 201210；3.中國科學(xué)院大學(xué)，北京 100049；4.中國科學(xué)院無錫高新微納傳感網(wǎng)工程技術(shù)研發(fā)中心，江蘇無錫 214135）

0 概述

基于傳感器陣列的波達(dá)方向（Direction of Arrival，DOA）估計被廣泛應(yīng)用于生產(chǎn)生活中［1］，例如雷達(dá)、聲吶、地震勘探、導(dǎo)航、聲源跟蹤等［2-4］?，F(xiàn)有DOA 估計算法大多假設(shè)被測信號是窄帶且噪聲服從高斯分布［5］，并利用多重信號分類［6-7］（Multiple Signal Classification，MUSIC）算法進(jìn)行估計，因此，在處理寬帶信號源時的性能受自然界聲波、地震波等因素［8-10］限制。

寬帶信號DOA 估計算法主要分為非相干信號子空間算法［11-12］和相干信號子空間算法［13-14］。非相干信號子空間算法僅分辨非相干信號，并且在低信噪比條件下估計性能不佳；相干信號子空間算法可用于相干信號，在低信噪比條件下估計準(zhǔn)確度高，但其聚焦過程運(yùn)算量較大且性能受波達(dá)角度預(yù)估計偏差的影響，波達(dá)角度預(yù)估計結(jié)果較差時會導(dǎo)致性能嚴(yán)重下降［15］。文獻(xiàn)［16］提出一種投影子空間正交性測試（TOPS）算法，利用多個頻點(diǎn)子空間的正交性實(shí)現(xiàn)寬帶信號DOA 估計，但其估計精度不高且易出現(xiàn)偽峰。文獻(xiàn)［17］提出修正的TOPS 算法，利用信號子空間投影有效剔除偽峰，但其在低信噪比條件下估計性能不佳。文獻(xiàn)［18］根據(jù)寬帶陣列導(dǎo)向矢量在法線方向上的頻率一致性，提出基于頻域時延補(bǔ)償?shù)腄OA 估計算法，該算法性能優(yōu)，但計算量較大。文獻(xiàn)［19］采用接收到的數(shù)據(jù)構(gòu)造聚焦矩陣以避免DOA 預(yù)估計，但估計性能依賴聚焦頻率的選取。文獻(xiàn)［20］利用陣型中理想的低秩Toeplitz 結(jié)構(gòu)，在快拍數(shù)不足的情況下實(shí)現(xiàn)DOA 估計，但該方法需要求解半定規(guī)劃問題且計算復(fù)雜度較高。文獻(xiàn)［21］提出信號子空間聚焦（Focusing Signal Subspace，F(xiàn)SS）算法，將參考頻率的信號子空間特征向量和其他頻率的信號子空間特征向量結(jié)合到Frobenius 范數(shù)的約束中，實(shí)現(xiàn)聚焦矩陣的構(gòu)造，該算法分辨率高且均方根誤差低，無需進(jìn)行初步的DOA 估計，但在短快拍條件下性能不佳。

本文提出一種改進(jìn)的信號子空間聚焦算法MFSS。根據(jù)子頻帶波長間隔與半波長的匹配度選取最佳參考頻率及子頻帶，減少聚焦過程的運(yùn)算量，同時將協(xié)方差矩陣平均處理為Hankel 矩陣并進(jìn)行奇異值分解重構(gòu)，降低噪聲及短快拍對協(xié)方差矩陣的影響。在此基礎(chǔ)上，利用信號子空間聚焦構(gòu)造最終的聚焦協(xié)方差矩陣，并通過Root-正交傳播算子實(shí)現(xiàn)DOA 估計。

1 寬帶信號模型

均勻線陣模型如圖1 所示。假設(shè)有P個相互獨(dú)立的遠(yuǎn)場寬帶相干信號以角度θ1，θ2，…，θP入射到M個陣元的均勻線陣上，陣元間隔d為信號最高頻率對應(yīng)的半波長。

圖1 均勻線陣模型Fig.1 Model of uniform linear array

第m個陣元接收的信號如式（1）所示：

其中：sp(t)為陣元接收的第p個信號；nm(t)為第m個陣元接收的噪聲；τmp為第p個信號到達(dá)第m個陣元較到達(dá)第1 個陣元的延遲。τmp如式（2）所示：

其中：c為信號傳播速度。

均勻線陣模型對接收數(shù)據(jù)xm(t)做離散傅里葉變換，將其劃分為J個子帶，快拍數(shù)為K，如式（3）所示：

矩陣形式的頻域陣列信號接收模型如式（4）所示：

其中：A(fj，θ)為陣列流型矩陣。列向量α(fj，θp)如式（5）所示：

X(fj)的協(xié)方差矩陣如式（6）所示：

其中：Rs(fj)、Rn(fj)分別為頻率上信號與噪聲的協(xié)方差矩陣。因此，Rs(fj)如式（7）所示：

2 基于MFSS 的寬帶信號DOA 估計

2.1 FSS 算法

FSS 算法的原理是結(jié)合參考頻率和其他頻率信號子空間的特征向量，在Frobenius 范數(shù)約束下構(gòu)造聚焦協(xié)方差矩陣，采用MUSIC 算法實(shí)現(xiàn)DOA 估計。

FSS 算法對R(fj)進(jìn)行特征分解，如式（8）所示：

其中：U(fj)為R(fj)的M×M維的特征向量矩陣；U(fj)=[e1(fj)，e2(fj)，…，eM(fj)]。由于信號與噪聲相互獨(dú)立，則信號子空間Us(fj)、噪聲子空間Un(fj)分別如式（9）、式（10）所示：

FSS 算法構(gòu)造非奇異矩陣T(fj)，將其作為聚焦矩陣，Us（fj）如式（11）所示：

為獲得最小的聚焦誤差，F(xiàn)SS 算法使用Frobenius范數(shù)約束信號子空間，如式（12）所示：

其中：T(fj)為Hermitian 矩陣。因此，聚焦矩陣T(fj)滿足式（13）：

將式（12）重寫為：

矩陣C如式（15）所示：

當(dāng)|Zii|=1，i=1，2，…，M時，式（17）可取最大值，此時：

則頻率fj處的聚焦協(xié)方差矩陣如式（20）所示：

最終的協(xié)方差矩陣RF是子頻帶的聚焦協(xié)方差矩陣的均值，如式（21）所示：

最后，通過MUSIC 算法實(shí)現(xiàn)DOA 估計。

2.2 MFSS 算法

MFSS 算法是將寬帶信號分為J段并進(jìn)行傅里葉變換，通過計算子頻帶波長間隔與半波長的插值選定參考頻率，并篩選出3 個子頻帶，將子頻帶的協(xié)方差矩陣處理為Hankel 矩陣，采用奇異值分解去噪并重構(gòu)協(xié)方差矩陣，利用信號子空間聚焦法構(gòu)造聚焦協(xié)方差矩陣并通過Root-正交傳播算子得到DOA估計值。MFSS 算法流程如圖2 所示。

圖2 MFSS 算法流程Fig.2 Procedure of MFSS algorithm

2.2.1 參考頻率及子頻帶選取

對于J個子頻帶，MFSS 算法計算頻率，其中c為速度，陣元間隔。在子帶1，2，…，J中，MFSS 算法搜索波長間隔最接近正整數(shù)j0=的子帶，其中mod 表示余數(shù)，fs為信號的采樣頻率。由于每個子頻帶信號的波長是不同的，因此MFSS算法需選取一個子頻帶與該波長間隔最為匹配。

MFSS 算法選取該子頻帶的中心頻率f0作為參考頻率，使用該子頻帶及臨近兩條子頻帶進(jìn)行后續(xù)的計算估計，臨近的兩條子頻帶頻率為fξ，ξ=1，2。

2.2.2 協(xié)方差矩陣重構(gòu)

在實(shí)際工程中，協(xié)方差矩陣由有限數(shù)量采樣條件下獲得接收數(shù)據(jù)的平均值構(gòu)成，且噪聲形式復(fù)雜。這些因素都會導(dǎo)致特征分解時信號子空間與噪聲子空間劃分模糊。針對該問題，本文采用Hankel 矩陣奇異值分解法對協(xié)方差矩陣進(jìn)行重構(gòu)。

本文對選取子頻帶的數(shù)據(jù)協(xié)方差矩陣R(fξ)次對角線及平行于次對角線直線上的元素進(jìn)行平均處理，如式（22）所示：

其中：k為對角線；k=0 為次對角線；0 ≤k

奇異值矩陣D中的大奇異值對應(yīng)信號分量，而小奇異值對應(yīng)噪聲分量，因此，保留奇異值矩陣D的前P個奇異值，后M-P個奇異值為0，以實(shí)現(xiàn)去噪的目的，處理后的矩陣為DP。奇異值矩陣?yán)肬、VT和DP構(gòu)建新的矩陣RP(fξ)，如式（25）所示：

2.2.3 信號子空間聚焦

本文對RP(fξ)及R(f0)特征進(jìn)行分解，求出fξ及f0處的信號子空間Us(fξ)及Us(f0)，對的乘積進(jìn)行奇異值分解，如式（26）所示：

2.2.4 Root-正交傳播算子

本文對協(xié)方差矩陣RMF進(jìn)行分塊，如式（28）所示：

其中：R1和R2分別為矩陣RMF的前P行及后M-P行，假設(shè)R1滿秩，則R2為R1的線性組合，故存在(M-P)×P維變換矩陣，如式（29）所示：

其中：IM-P為M-P維單位陣。Root-正交傳播算子對進(jìn)行正交化，如式（33）所示：

定義多項(xiàng)式為：

其中：z=ejw；p(z)=[1，z，…，zM-1]T。根據(jù)式（34）計算得到P個接近于單位元上的根，對于均勻線陣，單位元上的根可通過式（35）求解：

MFSS 算法的步驟主要分為：1）對陣列接收到的寬帶信號數(shù)據(jù)分段，并進(jìn)行離散傅里葉變換；2）選取波長間隔最接近子頻帶的中心頻率作為參考頻率，保留該子頻帶及其臨近兩條子頻帶；3）求得各頻點(diǎn)處的協(xié)方差矩陣R(fξ)，根據(jù)式（22）將其處理為Hankel 矩陣并利用奇異值分解重構(gòu)；4）對RP(fξ)及R(f0)進(jìn)行特征分解，根據(jù)各頻點(diǎn)處的信號子空間，構(gòu)造聚焦矩陣，計算聚焦協(xié)方差矩陣Rξ(f0)，與R(f0)進(jìn)行平均操作得到最終的協(xié)方差矩陣RMF；5）對協(xié)方差矩陣RMF分塊，根據(jù)式（32）求得噪聲子空間估計，計算式（35）的P個接近于單位元上的根，得到信號的DOA 估計。

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

為評估MFSS 算法的實(shí)用價值，本文分析MFSS、FSS、MTOPS［17］及LR-MUSIC 算法［20］的復(fù)雜度。MFSS算法的時間復(fù)雜度主要由以下5 項(xiàng)構(gòu)成：1）選取參考頻率及子頻帶時間復(fù)雜度O(J)；2）構(gòu)造協(xié)方差矩陣時間復(fù)雜度O(KM2)；3）重構(gòu)協(xié)方差矩陣RP時間復(fù)雜度O(M2+M3)；4）聚焦獲得最終協(xié)方差矩陣RMF時間復(fù)雜度O(M3)；5）構(gòu)造傳播算子多項(xiàng)式求根時間復(fù)雜度O(PM2)。因此，MFSS 算法總體時間復(fù)雜度約為O(J+PM2+M3)。FSS 算法需要J個子帶以實(shí)現(xiàn)聚焦協(xié)方差矩陣的構(gòu)造，并采用MUSIC 算法構(gòu)造空間譜，時間復(fù)雜度約為O(JKM2+JM3)。MTOPS算法構(gòu)造協(xié)方差矩陣并特征分解的時間復(fù)雜度約為O(JKM2+JM3)，計算判決矩陣的時間復(fù)雜度為O(2PM(M-P)+P2(M-P)(J-1))，求解矩陣跡時間復(fù)雜度為O(2P3(J-1))，MTOPS 算法總體時間復(fù)雜度較高。LR-MUSIC 算法構(gòu)造聚焦矩陣時間復(fù)雜度為O(JM3)，計算聚焦協(xié)方差矩陣時間復(fù)雜度為O(JKM2)，此時的時間復(fù)雜度已接近FSS 算法，在后續(xù)計算中最優(yōu)Toeplitz 矩陣通過半正定規(guī)劃求出，計算量大，且總體復(fù)雜度遠(yuǎn)大于MFSS 算法。與其他算法相比，MFSS 算法在時間復(fù)雜度上具有較大優(yōu)勢，更利于工程應(yīng)用。

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

本文對MFSS、FSS、MTOPS、LR-MUSIC 算法進(jìn)行仿真對比，以驗(yàn)證MFSS 算法在短快拍情況下的有效性。假設(shè)信號源數(shù)目已知，陣元數(shù)M=8，信號為頻率200～400Hz 的遠(yuǎn)場相干寬帶信號，陣元間距為信號中心頻率所對應(yīng)的半波長，噪聲為相互獨(dú)立且與信號無關(guān)的高斯白噪聲。DFT 點(diǎn)數(shù)為128，子頻帶數(shù)目J=20。LR-MUSIC算法是通過CVX工具包解決SDP問題。

4.1 實(shí)驗(yàn)1

假定2 個遠(yuǎn)場相干寬帶信號入射角分別為60°和80°，快拍數(shù)為50，本文進(jìn)行500 次蒙特卡洛實(shí)驗(yàn)。DOA 估計的均方根誤差（RRMSE）如式（36）所示：

其中：C為蒙特卡洛實(shí)驗(yàn)的數(shù)量。

在不同信噪比時4 種算法的均方根誤差對比如圖3 所示，其信噪比范圍從-20 dB 以間隔2 dB 升至5 dB。從圖3 可以看出，在短快拍條件下，4 種算法的均方根誤差均隨信噪比的增加而逐漸減小。在整個信噪比范圍內(nèi)，MFSS 算法的均方根誤差始終低于其他3 種算法。因此，在短快拍低信噪比條件下，MFSS 算法的估計誤差最小。

圖3 不同信噪比下4 種算法的均方根誤差對比Fig.3 Root mean square error comparison among four algorithms under different SNRs

圖4 不同信噪比下4 種算法估計成功率對比Fig.4 Estimation success rate comparison among four algorithms under different SNRs

4.2 實(shí)驗(yàn)2

本文考慮2 個入射角為θ1=60°、θ2=60°+Δθ的寬帶相干信號，信噪比為5 dB，快拍數(shù)為50。若，i=1，2，則算法可以成功分辨兩目標(biāo)。不同角度間隔下4種算法的估計成功率如圖5所示。

圖5 不同角度間隔下4 種算法的估計成功率對比Fig.5 Estimation success rate comparison among four algorithms under different angular separations

從圖5 可以看出，隨著寬帶信號源角度間隔增大，4 種算法的分辨性能明顯提高，當(dāng)兩信號源角度間隔為5°時，MFSS 和MTOPS 算法的估計成功率達(dá)到100%，而FSS 和LR-MUSIC 算法仍無法分辨兩信號源。在角度間隔小于5°時，MFSS 算法的估計成功率高于MTOPS 算法，MFSS 算法的成功估計角度間隔門限較FSS 算法降低了3°。因此，MFSS 算法在短快拍條件下能夠分辨角度間隔更小的信號源。

4.3 實(shí)驗(yàn)3

本文考慮2 個獨(dú)立寬帶相干信號源的入射角為θ1=60°、θ2=80°，信噪比為5 dB。在不同快拍數(shù)下4 種算法的均方根誤差對比如圖6 所示。從圖6 可以看出，快拍數(shù)從5 上升至50 時，MFSS 算法的均方根誤差最小。

圖6 不同快拍數(shù)下4 種算法的均方根誤差對比Fig.6 Root mean square error comparison among four algorithms under different snapshot values

在不同快拍數(shù)下4 種算法的估計成功率對比如圖7 所示。從圖7 可以看出，在快拍數(shù)大于40 時，MFSS、TOPS 和LR-MUSIC 算法的估計成功率較接近，逐漸趨于1。在短快拍條件下，MFSS 算法的估計成功率始終高于其他3 種算法。因此，在同條件下，MFSS 算法具有更高的估計精度。

圖7 不同快拍數(shù)下4 種算法的估計成功率對比Fig.7 Estimation success rate comparison among four algorithms under different snapshot values

4.4 實(shí)驗(yàn)4

本文分別考慮2 個獨(dú)立寬帶相干信號（入射角θ1=60°、θ2=80°）和3 個獨(dú)立寬帶相干信號（入射角θ1=40°、θ2=60°、θ3=80°）的情況，進(jìn)行500 次蒙特卡洛實(shí)驗(yàn)，得出算法的平均運(yùn)算時間。

不同算法的運(yùn)算時間對比如表1 所示。從表1可以看出，MFSS 算法的運(yùn)算時間最短，且遠(yuǎn)小于MTOPS 和LR-MUSIC 算法的運(yùn)算時間。相比FSS算法，在信源數(shù)為3 時，MFSS 算法的平均運(yùn)算時間降低了21.14%，更具實(shí)用性。

表1 不同算法的運(yùn)算時間對比Table 1 Computation time comparison among different algorithms s

5 結(jié)束語

本文提出一種無需角度預(yù)估計的信號子空間聚焦算法MFSS，利用奇異值的分布規(guī)律減少快拍數(shù)及噪聲對估計性能的影響，通過波長間隔與陣元間距的匹配度選取最佳參考頻點(diǎn)及子頻帶，降低運(yùn)算量。仿真結(jié)果表明，MFSS 算法在短快拍條件下能夠有效提高估計精度。后續(xù)將提高算法在復(fù)雜噪聲情況下的估計性能，使其適用于實(shí)際無線傳感網(wǎng)絡(luò)定位環(huán)境。