基于幾何序列分解與稀疏重構的DOA 估計

侯進，陳鑫強,3

（1.西南交通大學信息科學與技術學院智能感知智慧運維實驗室，四川 成都 611756；2.西南交通大學綜合交通大數(shù)據(jù)應用技術國家工程實驗室，四川 成都 611756；3.西南交通大學計算機與人工智能學院，四川 成都 611756）

0 引言

隨著通信、雷達等領域的發(fā)展，波達方向（DOA,direction of arrival）估計逐漸成為人們研究的熱點課題。DOA 估計的基本問題是利用陣列對空間中某一區(qū)域內感興趣的多個信號的空間位置進行估計，它面臨的挑戰(zhàn)主要集中在以下4 個方面：1) 在現(xiàn)實通信環(huán)境中，受高大建筑物和山丘的阻礙，電磁波在傳播時將發(fā)生反射，從而產(chǎn)生多徑效應，這將使部分DOA 估計算法失效；2) 出于成本等因素的考慮，實際陣列的陣元數(shù)往往小于目標信源數(shù)，因此，如何在欠定情況下進行DOA 估計成為一個難題；3) 在某些特定情況下，目標信源的來波方向可能相同，而目前許多DOA 估計算法無法對相同方向角的信源進行測向，因此，“角度兼并”問題成為DOA估計算法的另一挑戰(zhàn)；4) 多快拍數(shù)在一定程度上會降低DOA 估計算法的時效性，因此，如何在少快拍數(shù)下保證算法的高準確率和高分辨率，也成為DOA 估計算法待解決的一個問題。

目前，根據(jù)測向原理的不同，DOA 估計方法可分為3 類：基于特征結構子空間的方法、基于獨立成分分析（ICA,independent component analysis）的方法和基于稀疏重構的方法。最早的超分辨率DOA估計方法是著名的多信號分類（MUSIC,multiple signal classification）算法和基于旋轉不變技術的信號參數(shù)估計（ESPRIT,estimating signal parameter via rotational invariance technique）算法，它們都屬于特征結構子空間方法，由于相干信源會使空間協(xié)方差矩陣降秩，因此傳統(tǒng)的MUSIC 算法和ESPRIT 算法無法對相干信源進行測向。文獻[1-2]和文獻[3-4]分別提出了改進的MUSIC 算法和改進的ESPRIT算法用于對相干信號進行測向，但在欠定情況下該類算法性能將大幅度下降。受均勻線陣方向矢量結構的啟發(fā)，F(xiàn)riedlander 等[5]提出了一種將均勻圓陣轉換為虛擬均勻線陣的方法，并結合均勻線陣的解相干算法，如空間平滑算法[6-7]和矩陣重構算法[8-10]等，可以實現(xiàn)均勻圓陣對相干信號的測向，但在模式轉換的過程中，噪聲會進一步累積，因此該類算法在低信噪比（SNR,signal to noise ratio）下的性能較差。Ziskind 等[11]提出了一種基于確定性極大似然（DML,determined maximum likelihood）估計思想的DOA 估計方法，能夠實現(xiàn)均勻圓陣對相干信號的測向，但在求解過程中需要進行高維搜索，因此計算量較大，且不適合欠定情況下的DOA 估計。Zamani等[12]提出了一種改進的ICA 方法，該方法可以對混合信源進行盲分離，但在欠定情況和信源中存在相干信號時將失效?；谙∈柚貥嫷姆椒梢栽诘托旁氡群托旁粗写嬖谙喔尚盘柕那闆r下具有良好的性能，如Wang 等[13]和Trinh-Hoang 等[14]提出的正交匹配追蹤（OMP,orthogonal matching pursuit）算法及其改進算法、Soubies 等[15]和Zheng 等[16]提出的基于L2 范數(shù)的方法和基于稀疏貝葉斯的方法等，但在信源數(shù)較多和入射信源示向度間隔較小的情況下，該類算法的測向性能會有所下降。

針對現(xiàn)有算法的不足，本文提出了一種基于幾何序列分解和稀疏重構的方法用于DOA 估計。通過分析可以發(fā)現(xiàn)，陣列的等時間間隔采樣使接收序列成為多個復數(shù)等比序列的疊加，各個等比序列的公比蘊含了各個信源的頻率信息，而首項則蘊含了信源的幅值和相位信息。因此，本文首先利用幾何序列分解算法對混合序列中各個序列的公比和首項進行估計，估計的各個公比可用于信號重建，而估計的首項則可以構建各個信源的實際方向向量，即相干組標記。然后利用稀疏重構方法，即可實現(xiàn)對每個相干組進行DOA 估計。

1 系統(tǒng)模型

1.1 均勻圓陣模型

假設空間中均勻圓陣列的陣元數(shù)為M，半徑為r，陣元1 的空間位置為 [r,0]T，信號的來波方向相對于x軸的逆時針旋轉方向進行定義，如圖1 所示。

圖1 均勻圓陣示意

若從遠場發(fā)射而來的復數(shù)獨立信源數(shù)為k，則在t時刻，陣列的接收信號可表示為

其中，m=1,…,M，θi表示第i個獨立信源的波達角，λ表示信號波長。

1.2 相干信號模型

對于陣列接收的數(shù)據(jù)，可以利用相關系數(shù)來衡量接收信號之間的相關程度，其定義為

式(4)表明，當2 個信號相干時，它們之間只相差一個常復數(shù)。現(xiàn)假設第i個獨立信源中含有Pi個相干信號，則第i個獨立信源的接收數(shù)據(jù)可表示為

其中，a(θip) ∈CM、Bip和φip分別表示在第i個獨立信源中第p個相干信號的理論方向向量、幅值和相位。

令

此時，第i個獨立信源的接收數(shù)據(jù)可表示為

1.3 相干組標記

1.4 網(wǎng)格劃分

根據(jù)陣列相關理論可知，實際的陣列結構要求方向向量a(θ)與方向角一一對應，不能出現(xiàn)模糊現(xiàn)象，而入射信源的方向角相對于整個測向平面是稀疏的，因此，整個二維測向平面可以進行網(wǎng)格劃分，以構成方向角參數(shù)空間。在實際應用中，通常以1°為間隔進行網(wǎng)格劃分，改變方向角，使方向向量a(θ)在二維空間按方向角從小到大進行掃描，所形成的曲面稱為陣列流形，所構成的矩陣稱為過完備矩陣D∈CM×360，即

其中，Θ={1°,…,360°}表示方向角的參數(shù)空間。

2 基于幾何序列分解的實際方向向量估計

現(xiàn)考慮如下情況：假設均勻圓陣列的陣元數(shù)為M，入射到陣列的獨立信源數(shù)為k，第i個獨立信源中含有Pi個相干信號，信號接收機的采樣頻率為fs，快拍數(shù)為L，并假設采樣是從1 時刻開始，則陣列的接收數(shù)據(jù)可表示為

觀察式(13)可以發(fā)現(xiàn)，均勻圓陣第q個陣元的接收數(shù)據(jù)（觀測矩陣X的第q行數(shù)據(jù)）本質上是k個等比序列的疊加，其中，第i個等比序列的首項為，公比為，而首項和公比分別包含了第i個獨立信源的波達角信息和頻率信息。因此，若能估計出混合序列的k個公比，然后利用一個簡單的矩陣除法，就能將k個獨立信源的實際方向向量估計出來，從而用于DOA 估計。

Lee 等[17]提出了一種基于k維單純形轉換的幾何序列分解（GSD-ST,geometric sequence decomposition withk-simplexes transform）方法，可用于對疊加幾何序列的公比進行估計。k維單純形是指包含了k+1 個節(jié)點的凸多面體，例如，一維單純形就是包含了2 個節(jié)點的線段，二維單純形就是三角形，三維單純形就是四面體，以此類推。利用GSD-ST方法估計公比的思路主要分為4 個步驟：構造搜索空間；構造k維單純形序列；構造聯(lián)合單純形；構造多項式方程。而多項式方程的根即待求的各個序列公比的估計。

在構造k維單純形之前，需要對原始的疊加序列按照一定方式進行采樣，以構成k維單純形的頂點集合，也稱為搜索空間，其構造規(guī)則按以下方式進行。

給定自然數(shù)ic、k和φk，以φk為起始坐標、i Ick為索引變換按字典順序對混合序列xq進行連續(xù)采樣，所構成的集合稱為k維搜索空間，其中，I k表示k維的全1 向量，xq表示第q個陣元接收的混合序列。例如，若ic=1，k=2，φk= [1,5]，則所構成的k維搜索空間為

由k維單純形的定義可知，一個k維單純形具有k+1 個頂點，現(xiàn)假設其中一個頂點為坐標原點，剩余k個頂點由搜索空間中的k個連續(xù)頂點構成，其定義為

需要指出的是，在構造本文要求的k維單純形時，從搜索空間中選取的k個頂點必須是連續(xù)的，且在構造k維單純形序列時，坐標j的取值也需要滿足連續(xù)遞增的條件。

接著，從構造的k維單純形序列中挑選出任意2 個連續(xù)的單純形，將其坐標按順序排列，以構成一個具有k+2 個頂點的聯(lián)合單純形C∈Ck×(k+2)，定義為

從式(17)可以看出，聯(lián)合單純形C的k+1 個非原點坐標實際上由搜索空間的k+1 個連續(xù)頂點坐標構成；然后，對這k+1 個連續(xù)頂點按組合的順序進行采樣，以構成一個含有k+1 個新的k維單純形的序列。例如，若ic= 1,k= 2,φk= [1,5],j=1，則所構成的聯(lián)合單純形C和含有k+1 個新的k維單純形的序列為

最后，需要構造一個k階多項式，并且這個k階多項式的k+1 個系數(shù)由序列中單純形的體積通過運算得到。其中，k維單純形的體積定義為

將k維單純形序列的體積運算規(guī)則定義為

其中，表示公比的估計。

Lee 等[17]證明，給定混合序列xq和獨立信源數(shù)k，當采樣間隔ic=1時，對于任意的自然數(shù)φk和j，k階多項式的系數(shù)集合是唯一的，并且由這個系數(shù)集合構成的k階多項式的k個根即混合序列中k個公比的估計。由于系數(shù)集合在上述采樣方式下與常數(shù)φk和j無關，因此，可以把系數(shù)集合記為v(k,xq)。需要注意的是，為了使構造的多項式是人們期望的結構，也就是要保證這個多項式的唯一性，在構造搜索空間時，一定要以1×Ik為索引變換對混合序列進行采樣，即采樣間隔一定要等于1，而采樣起始坐標φk和聯(lián)合多面體C的索引值j可以任取。同時，Lee 等[17]還證明，給定混合序列xq和獨立信源數(shù)k，獲取的最小樣本集合由混合序列xq的2k個連續(xù)采樣點構成，即在估計k個獨立信源的實際方向向量時，所需要的最小快拍數(shù)L為2k，這也使該算法可以在極少快拍數(shù)下對DOA 進行有效估計。

通過以上分析可以發(fā)現(xiàn)，在對混合序列中各個序列的公比進行估計時，需要預先知道混合序列中的序列個數(shù)，也就是獨立信源的個數(shù)k。Lee等[17]還提出了一種利用單純形對獨立信源數(shù)k進行估計的方法，但該方法需要進行多維搜索，計算量較大，因此會大幅度降低DOA 估計的時效性。Bazzi 等[18]提出了一種改進的最小描述長度（MDL,minimum description length）準則，可以實現(xiàn)在信噪比較低和快拍數(shù)較少的情況下對獨立信源數(shù)k進行有效估計，且算法時效性較好，因此，可以使用該方法對k進行估計。

受現(xiàn)實擾動的影響，陣列在接收數(shù)據(jù)時會產(chǎn)生一定的噪聲，因此，為了提升DOA 估計的性能，需要對接收數(shù)據(jù)進行去噪。Hansen[19]提出了一種基于k截斷式奇異值分解（SVD,singular value decomposition）的去噪方法，該方法的主要思想是通過迭代的方法，不斷提取原序列的k個主特征值和它們對應的特征向量，直到重構的序列收斂，從而達到濾除噪聲的目的。該方法的去噪過程主要分為以下5 個步驟，這里假設去噪前和去噪后的序列分別為xq和，序列的索引值從1 開始，且快拍數(shù)為L。

其中，矩陣U和矩陣V分別表示左奇異矩陣和右奇異矩陣，矩陣S為對角矩陣。

步驟3利用矩陣Q的k個主特征值和對應的特征向量重構矩陣Q

其中，U[:,1:k]表示由矩陣U的所有行和前k列構成的矩陣。

步驟4通過平均值的方法，將矩陣Q轉換為序列

其中，當l≤m時，z=l，否則，z=L-l+1。

步驟5重復以上4 個步驟，直到序列xq收斂到序列，收斂后的序列即濾除噪聲后的序列。

以上就是估計k個實際方向向量的整體思路，現(xiàn)將算法流程總結如下。

算法1基于去噪GSD-ST 的實際方向向量估計算法

3 基于稀疏重構的DOA 估計

根據(jù)上述分析可知，通過網(wǎng)格劃分構造了一個過完備矩陣D，并利用過完備矩陣中的稀疏基對第i個相干組的實際方向向量進行了稀疏表示，將第i個相干組中多徑信號的DOA 估計問題轉換為空間譜Vi的稀疏重構問題，若對每個相干組都進行稀疏重構，得到k張空間譜(V1,…,Vk)，即可估計出所有入射信源的波達角信息。

在求解第i個相干組的稀疏矩陣Vi時，由于稀疏度Pi事先未知，且實際方向向量的維度遠小于矩陣Vi的維度，因此式(26)具有無數(shù)個解。根據(jù)壓縮理論知識可知，如果過完備矩陣D滿足約束等距性（RIP,restricted isometry property）條件，即可利用矩陣Vi的稀疏性，將式(26)轉換為一個組合優(yōu)化問題進行求解?，F(xiàn)對該條件進行簡單驗證，由于

在求解式(28)時，可以通過正則化對這個非線性凸優(yōu)化問題進行求解，其代價函數(shù)為

其中，λ為正則化系數(shù)。

在實際仿真時，為了高效地找到全局最優(yōu)解，可以利用二階錐規(guī)劃（SOCP,second order cone programming）算法在內點法的框架下對式(29)中的優(yōu)化表達式進行求解，基于此，可以通過引入輔助向量的方式，將式(29)轉換為以下優(yōu)化表達式進行求解

需要注意的是，正則化系數(shù)λ的選取對稀疏重構的效果有著舉足輕重的作用，文獻[14]給出了正則化系數(shù)λ的選取規(guī)則，由于本文通過k截斷式SVD 算法對原觀測矩陣進行了去噪處理，因此系數(shù)λ的值不應取得太大，在實際仿真時，系數(shù)λ在區(qū)間[0.1,3]內進行取值都能達到很好的效果。

以上就是本文所提算法的全部分析過程，現(xiàn)將算法流程總結如下。

算法2基于去噪GSD-ST 和稀疏重構的DOA估計算法

1) 利用算法1 得到獨立信源數(shù)k和實際方向矩陣的估計

2) 以1°為間隔對空間進行網(wǎng)格劃分，構造過完備矩陣D；

5) 將稀疏矩陣Vi的元素值畫成空間譜，通過譜峰搜索返回第i個相干組的DOA 估計值；

6)i=i+1，如果i≤k，返回步驟4)，否則算法結束。

4 性能分析

由于本文使用的信源數(shù)估計算法為改進的MDL 算法，因此，基于GSD-ST 的實際方向向量估計算法最多能夠估計出M個獨立信源的實際方向向量；并且，一般情況下入射信源的幅值滿足隨機分布，對于每一個相干組，稀疏重構算法能夠估計出M–1 個相干信號[20]，因此，本文所提算法最多能夠估計出M(M–1)個入射信源，即能夠實現(xiàn)在欠定情況下對相干信號的DOA 估計；在第2 節(jié)的分析中，本文對相干組進行了拆分，即每個相干組的DOA 估計互不影響，因此，本文所提算法能夠解決“角度兼并”問題；此外，在估計實際方向向量時，所需要的最小快拍數(shù)L=2k，因此，本文所提算法可以在極少數(shù)快拍下實現(xiàn)DOA 估計。

5 仿真分析

為了驗證本文所提算法的性能，本節(jié)將其與非線性快速獨立成分分析（FastICA,fast independent component analysis）算法[20]、OMP 算法[13]、MUSIC算法[1]和DML 算法[11]進行對比，并通過5 個實驗進行說明。

實驗1相干組拆分性能分析

實驗條件：均勻圓陣的陣元數(shù)M=3，半徑r=0.58 m，信號的中心頻率f=500 MHz，相干組個數(shù)k=3，正確示向度分別為[32°]、[210°]、[89°]，信噪比SNR=[-18,18] dB，快拍數(shù)L=200，蒙特卡羅實驗200 次。相干組拆分效果通過改進的行元素優(yōu)勢指標[21]進行分析，定義為

圖2 為非線性FastICA 算法與本文所提算法的盲分離效果對比。從圖2 中可以看出，相比于非線性FastICA 算法，本文所提算法在低信噪比下的p指標系數(shù)都小于2，這表明本文所提算法可以在低信噪比下對相干組進行有效拆分，且拆分效果優(yōu)于非線性FastICA 算法。

圖2 2 種算法盲分離效果對比

為了進一步說明本文所提算法拆分相干組的效果，圖3 和圖4 給出了當信噪比為-18 dB 時相干組拆分前后的數(shù)據(jù)頻譜曲線。

圖3 拆分前的數(shù)據(jù)頻譜曲線

圖4 拆分后的數(shù)據(jù)頻譜曲線

對比圖3 和圖4 可以看出，混合數(shù)據(jù)的頻譜中有許多偽峰，而本文所提算法對混合數(shù)據(jù)進行盲分離后，其頻譜曲線濾除了由噪聲引起的偽峰，并且能夠準確恢復出信號的頻率信息，這也驗證了算法1的正確性。

實驗2 欠定情況下對相干信號的DOA 估計

實驗條件：均勻圓陣的陣元數(shù)M=3，半徑r=0.58 m，信號的中心頻率f=500 MHz，相干組個數(shù)k=3，相干系數(shù)分別為=-0 .416+0.909i、= 0.921+0.389i 和= 0.853 - 0.522i，正確示向度分別為[32°,328°]、[30°,210°]、[89°,272°]，信噪比SNR=10 dB，快拍數(shù)L=200。

在上述實驗條件下，改進的MDL 算法估計出的獨立信源（相干組）數(shù)k=3，因此，MUSIC算法和DML 算法只能估計出3 個信號的示向度，為了驗證本文所提算法在欠定情況下對相干信號的DOA 估計性能，將其與OMP 算法進行了對比。圖5 和表1 分別為OMP 算法與本文所提算法在欠定情況下對相干信號的測向空間譜和對應的測向結果。

圖5 2 種算法在欠定情況下對相干信號的測向空間譜

表1 2 種算法在欠定情況下對相干信號的測向結果

從圖5 和表1 可以看出，在欠定情況下，OMP算法估計出的DOA 結果相較正確結果的誤差最大達到了3°，并且由于30°和32°這2 個示向度較接近，OMP 算法不能有效分辨，從而出現(xiàn)了3°這個錯誤估計結果；本文所提算法不僅能夠有效拆分出3 個相干組，且DOA 估計結果全部正確，這也驗證了本文所提算法所能準確估計的最大信源數(shù)為M(M–1)的結論。

實驗3極少快拍數(shù)下的測向性能

實驗條件均勻圓陣的陣元數(shù)M=9，半徑r=0.58 m，信號的中心頻率f=500 MHz，獨立信源數(shù)k=3，正確示向度分別為[15°]、[167°]和[251°]，信噪比SNR=30 dB，快拍數(shù)L=6。

在極少快拍數(shù)的情況下，本節(jié)對比了OMP 算法、MUSIC 算法和DML 算法，由于DML 算法需要進行高維搜索，因此沒有畫出DML 算法的空間譜，只記錄了對應的測向結果。圖6 和表2 分別為對比算法與本文所提算法在極少快拍數(shù)下的測向空間譜和對應的測向結果。

圖6 3 種算法在極少快拍數(shù)下的測向空間譜

表2 4 種算法在極少快拍數(shù)下的測向結果

從圖6 和表2 可以看出，在極少快拍數(shù)的情況下，OMP 算法、MUSIC 算法和DML 算法在測向時會有2°～3°的測向誤差，而本文所提算法能夠精確地對入射信號進行DOA 估計。在實驗過程中還發(fā)現(xiàn)，快拍數(shù)的減小會使噪聲對測向的影響變得很大，因此，極少快拍數(shù)測向只適用于信噪比較高的情況，雖然現(xiàn)實測向情況較復雜，但本文所提算法為極少快拍數(shù)測向提供了一定的理論和仿真依據(jù)。

實驗4入射信源角度差對測向成功率的影響

實驗條件均勻圓陣的陣元數(shù)M=9，半徑r=0.58 m，信號的中心頻率f=500 MHz，獨立信源數(shù)k=2，正確示向度分別為 [123°] 和[123°+Δθ(Δθ=0°,…,6°)]，信噪比SNR=10 dB，快拍數(shù)L=200，蒙特卡羅實驗200 次。

圖7 為4 種算法測向成功率隨信號角度差變化的曲線。從圖7 可以看出，OMP 算法、MUSIC 算法和DML 算法都無法分辨出角度差為0°的2 個入射信號，并且OMP 算法和MUSIC 算法在入射信源夾角超過4°時才能達到很好的分辨效果。本文所提算法由于對相干組進行了拆分，相干組之間的測向互不影響，因此可以對角度差較小的入射信號進行精確分辨。

圖7 4 種算法測向成功率隨信號角度差變化的曲線

實驗5信噪比對測向成功率和測向精度的影響

實驗條件均勻圓陣的陣元數(shù)M=9，半徑r=0.58 m，信號的中心頻率f=500 MHz，獨立信源數(shù)k=3，正確示向度分別為[52°]、[189°]和[297°]，快拍數(shù)L=200，蒙特卡羅實驗200 次。其中，測向精度采用均方誤差根（RMSE,root mean square error）進行定義

其中，Q表示測向正確示向度的次數(shù)，即測向誤差結果在3°范圍內的個數(shù)；表示在第q次測向時對第i個示向度的估計值。

圖8 和圖9 分別給出了4 種測向算法的成功率和RMSE 隨信噪比變化的曲線。從圖8 可以看出，當信噪比大于-16 dB 時，本文所提算法的測向成功率可達90%以上，且高于OMP 算法和MUSIC 算法；當信噪比大于-14 dB 時，本文所提算法的測向成功率可達99%以上，這充分說明了本文所提算法在低信噪比下能夠對多個信號進行有效的DOA 估計。從圖9 可以看出，當信噪比大于-16 dB 時，本文所提算法的測向精度達到了0.5°以上，且對應的精度曲線收斂到了0°，而OMP算法、MUSIC 算法和DML 算法對應的測向精度曲線分別收斂到了0.8°、1.2°和0.7°。綜合來看，相比于其他3 種測向算法，本文所提算法的抗噪能力和測向精度更優(yōu)。

圖8 4 種算法測向成功率隨信噪比變化的曲線

圖9 4 種算法測向精度隨信噪比變化的曲線

6 結束語

本文以提升均勻圓陣對相干信號的測向性能為目標，提出了一種去噪GSD-ST 與稀疏重構相結合的DOA 估計算法，去噪GSD-ST 算法保證了對相干組拆分和估計實際方向矩陣的有效性和穩(wěn)定性，稀疏重構算法則利用估計的實際方向矩陣對每個相干組進行DOA 估計。仿真實驗表明，在欠定情況和入射信源存在相干信號時，本文所提算法的測向準確率和抗擾能力都更優(yōu)。當入射信源角度差較小和信源數(shù)較多時，本文所提算法的測向優(yōu)勢更明顯；同時，本文所提算法為極少快拍數(shù)測向提供了理論和仿真依據(jù)。接下來，筆者將進一步探索如何在低信噪比下利用極少快拍數(shù)進行測向的問題，以進一步提高算法的測向性能。