賈思宇,路 茗,丁華澤,陳 明,趙魯陽
(1.中國科學院上海微系統與信息技術研究所無線傳感網與通信重點實驗室,上海 200050;2.上??萍即髮W信息科學與技術學院,上海 201210;3.中國科學院大學,北京 100049;4.中國科學院無錫高新微納傳感網工程技術研發中心,江蘇無錫 214135)
基于傳感器陣列的波達方向(Direction of Arrival,DOA)估計被廣泛應用于生產生活中[1],例如雷達、聲吶、地震勘探、導航、聲源跟蹤等[2-4]?,F有DOA 估計算法大多假設被測信號是窄帶且噪聲服從高斯分布[5],并利用多重信號分類[6-7](Multiple Signal Classification,MUSIC)算法進行估計,因此,在處理寬帶信號源時的性能受自然界聲波、地震波等因素[8-10]限制。
寬帶信號DOA 估計算法主要分為非相干信號子空間算法[11-12]和相干信號子空間算法[13-14]。非相干信號子空間算法僅分辨非相干信號,并且在低信噪比條件下估計性能不佳;相干信號子空間算法可用于相干信號,在低信噪比條件下估計準確度高,但其聚焦過程運算量較大且性能受波達角度預估計偏差的影響,波達角度預估計結果較差時會導致性能嚴重下降[15]。文獻[16]提出一種投影子空間正交性測試(TOPS)算法,利用多個頻點子空間的正交性實現寬帶信號DOA 估計,但其估計精度不高且易出現偽峰。文獻[17]提出修正的TOPS 算法,利用信號子空間投影有效剔除偽峰,但其在低信噪比條件下估計性能不佳。文獻[18]根據寬帶陣列導向矢量在法線方向上的頻率一致性,提出基于頻域時延補償的DOA 估計算法,該算法性能優,但計算量較大。文獻[19]采用接收到的數據構造聚焦矩陣以避免DOA 預估計,但估計性能依賴聚焦頻率的選取。文獻[20]利用陣型中理想的低秩Toeplitz 結構,在快拍數不足的情況下實現DOA 估計,但該方法需要求解半定規劃問題且計算復雜度較高。文獻[21]提出信號子空間聚焦(Focusing Signal Subspace,FSS)算法,將參考頻率的信號子空間特征向量和其他頻率的信號子空間特征向量結合到Frobenius 范數的約束中,實現聚焦矩陣的構造,該算法分辨率高且均方根誤差低,無需進行初步的DOA 估計,但在短快拍條件下性能不佳。
本文提出一種改進的信號子空間聚焦算法MFSS。根據子頻帶波長間隔與半波長的匹配度選取最佳參考頻率及子頻帶,減少聚焦過程的運算量,同時將協方差矩陣平均處理為Hankel 矩陣并進行奇異值分解重構,降低噪聲及短快拍對協方差矩陣的影響。在此基礎上,利用信號子空間聚焦構造最終的聚焦協方差矩陣,并通過Root-正交傳播算子實現DOA 估計。
均勻線陣模型如圖1 所示。假設有P個相互獨立的遠場寬帶相干信號以角度θ1,θ2,…,θP入射到M個陣元的均勻線陣上,陣元間隔d為信號最高頻率對應的半波長。

圖1 均勻線陣模型Fig.1 Model of uniform linear array
第m個陣元接收的信號如式(1)所示:

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

其中:c為信號傳播速度。
均勻線陣模型對接收數據xm(t)做離散傅里葉變換,將其劃分為J個子帶,快拍數為K,如式(3)所示:

矩陣形式的頻域陣列信號接收模型如式(4)所示:

其中:A(fj,θ)為陣列流型矩陣。列向量α(fj,θp)如式(5)所示:

X(fj)的協方差矩陣如式(6)所示:

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

FSS 算法的原理是結合參考頻率和其他頻率信號子空間的特征向量,在Frobenius 范數約束下構造聚焦協方差矩陣,采用MUSIC 算法實現DOA 估計。
FSS 算法對R(fj)進行特征分解,如式(8)所示:

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

FSS 算法構造非奇異矩陣T(fj),將其作為聚焦矩陣,Us(fj)如式(11)所示:

為獲得最小的聚焦誤差,FSS 算法使用Frobenius范數約束信號子空間,如式(12)所示:

其中:T(fj)為Hermitian 矩陣。因此,聚焦矩陣T(fj)滿足式(13):

將式(12)重寫為:

矩陣C如式(15)所示:



當|Zii|=1,i=1,2,…,M時,式(17)可取最大值,此時:

則頻率fj處的聚焦協方差矩陣如式(20)所示:

最終的協方差矩陣RF是子頻帶的聚焦協方差矩陣的均值,如式(21)所示:

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

圖2 MFSS 算法流程Fig.2 Procedure of MFSS algorithm
2.2.1 參考頻率及子頻帶選取
對于J個子頻帶,MFSS 算法計算頻率,其中c為速度,陣元間隔。在子帶1,2,…,J中,MFSS 算法搜索波長間隔最接近正整數j0=的子帶,其中mod 表示余數,fs為信號的采樣頻率。由于每個子頻帶信號的波長是不同的,因此MFSS算法需選取一個子頻帶與該波長間隔最為匹配。
MFSS 算法選取該子頻帶的中心頻率f0作為參考頻率,使用該子頻帶及臨近兩條子頻帶進行后續的計算估計,臨近的兩條子頻帶頻率為fξ,ξ=1,2。
2.2.2 協方差矩陣重構
在實際工程中,協方差矩陣由有限數量采樣條件下獲得接收數據的平均值構成,且噪聲形式復雜。這些因素都會導致特征分解時信號子空間與噪聲子空間劃分模糊。針對該問題,本文采用Hankel 矩陣奇異值分解法對協方差矩陣進行重構。
本文對選取子頻帶的數據協方差矩陣R(fξ)次對角線及平行于次對角線直線上的元素進行平均處理,如式(22)所示:

其中:k為對角線;k=0 為次對角線;0 ≤k 奇異值矩陣D中的大奇異值對應信號分量,而小奇異值對應噪聲分量,因此,保留奇異值矩陣D的前P個奇異值,后M-P個奇異值為0,以實現去噪的目的,處理后的矩陣為DP。奇異值矩陣利用U、VT和DP構建新的矩陣RP(fξ),如式(25)所示: 2.2.3 信號子空間聚焦 本文對RP(fξ)及R(f0)特征進行分解,求出fξ及f0處的信號子空間Us(fξ)及Us(f0),對的乘積進行奇異值分解,如式(26)所示: 2.2.4 Root-正交傳播算子 本文對協方差矩陣RMF進行分塊,如式(28)所示: 其中:R1和R2分別為矩陣RMF的前P行及后M-P行,假設R1滿秩,則R2為R1的線性組合,故存在(M-P)×P維變換矩陣,如式(29)所示: 其中:IM-P為M-P維單位陣。Root-正交傳播算子對進行正交化,如式(33)所示: 定義多項式為: 其中:z=ejw;p(z)=[1,z,…,zM-1]T。根據式(34)計算得到P個接近于單位元上的根,對于均勻線陣,單位元上的根可通過式(35)求解: MFSS 算法的步驟主要分為:1)對陣列接收到的寬帶信號數據分段,并進行離散傅里葉變換;2)選取波長間隔最接近子頻帶的中心頻率作為參考頻率,保留該子頻帶及其臨近兩條子頻帶;3)求得各頻點處的協方差矩陣R(fξ),根據式(22)將其處理為Hankel 矩陣并利用奇異值分解重構;4)對RP(fξ)及R(f0)進行特征分解,根據各頻點處的信號子空間,構造聚焦矩陣,計算聚焦協方差矩陣Rξ(f0),與R(f0)進行平均操作得到最終的協方差矩陣RMF;5)對協方差矩陣RMF分塊,根據式(32)求得噪聲子空間估計,計算式(35)的P個接近于單位元上的根,得到信號的DOA 估計。 為評估MFSS 算法的實用價值,本文分析MFSS、FSS、MTOPS[17]及LR-MUSIC 算法[20]的復雜度。MFSS算法的時間復雜度主要由以下5 項構成:1)選取參考頻率及子頻帶時間復雜度O(J);2)構造協方差矩陣時間復雜度O(KM2);3)重構協方差矩陣RP時間復雜度O(M2+M3);4)聚焦獲得最終協方差矩陣RMF時間復雜度O(M3);5)構造傳播算子多項式求根時間復雜度O(PM2)。因此,MFSS 算法總體時間復雜度約為O(J+PM2+M3)。FSS 算法需要J個子帶以實現聚焦協方差矩陣的構造,并采用MUSIC 算法構造空間譜,時間復雜度約為O(JKM2+JM3)。MTOPS算法構造協方差矩陣并特征分解的時間復雜度約為O(JKM2+JM3),計算判決矩陣的時間復雜度為O(2PM(M-P)+P2(M-P)(J-1)),求解矩陣跡時間復雜度為O(2P3(J-1)),MTOPS 算法總體時間復雜度較高。LR-MUSIC 算法構造聚焦矩陣時間復雜度為O(JM3),計算聚焦協方差矩陣時間復雜度為O(JKM2),此時的時間復雜度已接近FSS 算法,在后續計算中最優Toeplitz 矩陣通過半正定規劃求出,計算量大,且總體復雜度遠大于MFSS 算法。與其他算法相比,MFSS 算法在時間復雜度上具有較大優勢,更利于工程應用。 本文對MFSS、FSS、MTOPS、LR-MUSIC 算法進行仿真對比,以驗證MFSS 算法在短快拍情況下的有效性。假設信號源數目已知,陣元數M=8,信號為頻率200~400Hz 的遠場相干寬帶信號,陣元間距為信號中心頻率所對應的半波長,噪聲為相互獨立且與信號無關的高斯白噪聲。DFT 點數為128,子頻帶數目J=20。LR-MUSIC算法是通過CVX工具包解決SDP問題。 假定2 個遠場相干寬帶信號入射角分別為60°和80°,快拍數為50,本文進行500 次蒙特卡洛實驗。DOA 估計的均方根誤差(RRMSE)如式(36)所示: 其中:C為蒙特卡洛實驗的數量。 在不同信噪比時4 種算法的均方根誤差對比如圖3 所示,其信噪比范圍從-20 dB 以間隔2 dB 升至5 dB。從圖3 可以看出,在短快拍條件下,4 種算法的均方根誤差均隨信噪比的增加而逐漸減小。在整個信噪比范圍內,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 本文考慮2 個入射角為θ1=60°、θ2=60°+Δθ的寬帶相干信號,信噪比為5 dB,快拍數為50。若,i=1,2,則算法可以成功分辨兩目標。不同角度間隔下4種算法的估計成功率如圖5所示。 圖5 不同角度間隔下4 種算法的估計成功率對比Fig.5 Estimation success rate comparison among four algorithms under different angular separations 從圖5 可以看出,隨著寬帶信號源角度間隔增大,4 種算法的分辨性能明顯提高,當兩信號源角度間隔為5°時,MFSS 和MTOPS 算法的估計成功率達到100%,而FSS 和LR-MUSIC 算法仍無法分辨兩信號源。在角度間隔小于5°時,MFSS 算法的估計成功率高于MTOPS 算法,MFSS 算法的成功估計角度間隔門限較FSS 算法降低了3°。因此,MFSS 算法在短快拍條件下能夠分辨角度間隔更小的信號源。 本文考慮2 個獨立寬帶相干信號源的入射角為θ1=60°、θ2=80°,信噪比為5 dB。在不同快拍數下4 種算法的均方根誤差對比如圖6 所示。從圖6 可以看出,快拍數從5 上升至50 時,MFSS 算法的均方根誤差最小。 圖6 不同快拍數下4 種算法的均方根誤差對比Fig.6 Root mean square error comparison among four algorithms under different snapshot values 在不同快拍數下4 種算法的估計成功率對比如圖7 所示。從圖7 可以看出,在快拍數大于40 時,MFSS、TOPS 和LR-MUSIC 算法的估計成功率較接近,逐漸趨于1。在短快拍條件下,MFSS 算法的估計成功率始終高于其他3 種算法。因此,在同條件下,MFSS 算法具有更高的估計精度。 圖7 不同快拍數下4 種算法的估計成功率對比Fig.7 Estimation success rate comparison among four algorithms under different snapshot values 本文分別考慮2 個獨立寬帶相干信號(入射角θ1=60°、θ2=80°)和3 個獨立寬帶相干信號(入射角θ1=40°、θ2=60°、θ3=80°)的情況,進行500 次蒙特卡洛實驗,得出算法的平均運算時間。 不同算法的運算時間對比如表1 所示。從表1可以看出,MFSS 算法的運算時間最短,且遠小于MTOPS 和LR-MUSIC 算法的運算時間。相比FSS算法,在信源數為3 時,MFSS 算法的平均運算時間降低了21.14%,更具實用性。 表1 不同算法的運算時間對比Table 1 Computation time comparison among different algorithms s 本文提出一種無需角度預估計的信號子空間聚焦算法MFSS,利用奇異值的分布規律減少快拍數及噪聲對估計性能的影響,通過波長間隔與陣元間距的匹配度選取最佳參考頻點及子頻帶,降低運算量。仿真結果表明,MFSS 算法在短快拍條件下能夠有效提高估計精度。后續將提高算法在復雜噪聲情況下的估計性能,使其適用于實際無線傳感網絡定位環境。











3 算法復雜度分析
4 仿真實驗與分析
4.1 實驗1



4.2 實驗2

4.3 實驗3


4.4 實驗4

5 結束語