基于改進典范分解的嵌套陣列DOA估計

程思備，駱驍，蔣博籌，王玉婷，吳寰宇

（重慶信息通信研究院，重慶 401336）

1 引言

在陣列信號處理中，更高的自由度（degree of freedoms，DOF）往往能夠帶來更高的輸出信噪比、更低的旁瓣、更低的克拉美羅下界（Cramer-Rao lower bound，CRLB）等。與均勻線陣（uniformlinear array，ULA）相比，非均勻線陣（nonuniform linear array，NLA）因為其稀疏特性能夠在使用相同陣元數的情況下，獲得更高的 DOF[1]。常見的非均勻線陣有最小不冗余陣（minimum redundancy array，MRA）[2]、互質陣列（co-prime array）[3-4]以及嵌套陣列（nested array）[5]。

參考文獻[2]結合 MRA結構提出了基于擴展矩陣的算法，該算法的確能達到增加DOF的目的，從而獲得更優異的處理性能。但是在給定陣元數的情況下，MRA的布陣方式沒有相應的解析表達式，只能通過迭代算法得到其陣元位置，因此該陣列通用性較差。2010年，嵌套陣列的概念和數學模型首次被提出，該陣列由兩個或者多個陣元間隔不同的均勻線陣連接而成，當陣元數和階數確定后其陣元位置表達式可由固定計算式直接給出[5]，近幾年已經成為非均勻陣列領域研究的熱點。參考文獻[6]提出了超級嵌套陣列（super nested array）的布陣方式，該方式保留了嵌套陣列增加自由度優點的同時也能夠減少陣元間的互耦。對于嵌套陣列的DOA估計，現有參考文獻采用多重信號分類 （multiple signal classification，MUSIC）[7]、傳播因子與 MUSIC 結合[8]、酉變換矩陣重構與MUSIC結合[9]等算法，上述算法需要在角度域進行峰值搜索，計算量較大。傳統的旋轉平移不變技術（estimating signal parameters via rotational invariance technique，ESPRIT）不需要進行峰值搜索，它利用均勻陣列的導向量矩陣為范德蒙德矩陣（Vandermonde matrix），其具有全局平移不變的特性，將陣列分為去頂和去底的子陣列，而通過求解子陣列協方差矩陣特征向量之間的旋轉因子來完成 DOA估計[10]，但是全局平移不變特性不具有普遍性，只有均勻陣列具有該特性，大多數陣列只有局部平移不變特性，因此傳統的 ESPRIT算法不適用于嵌套陣列的DOA估計。

參考文獻[11-13]利用了非均勻線陣局部平移不變特性（將非均勻線陣分為若干個的子陣列，每個子陣列為均勻線陣，其導向量矩陣均為范德摩爾矩陣，具有平移不變性），將子陣列導向量矩陣構造轉化為三階張量的構造，最后利用典范分解求解三階張量的水平切片，得到非均勻線陣的導向量矩陣，最后利用最小二范數法則完成DOA估計。ESPRIT算法與基于典范分解的DOA估計算法類似，都利用范德摩爾矩陣平移不變特性，不同的是ESPRIT算法利用全局特性而DOA估計算法利用的是局部特性，因此基于典范分解的DOA估計方法被認為是ESPRIT算法的延伸，被定義為“Multiple Invariance ESPRIT”[13]。但是在參考文獻[13]中，快拍（Snapshot）信號模型不含噪聲成分，作者指出該算法只適用于無噪環境的假設前提，而在實際環境中噪聲普遍存在，通過仿真發現在有背景噪聲的條件下，算法的魯棒性急劇下降，甚至出現無解。如果將此算法應用于實際工程，必須對其進行改進優化。

本文結合含有噪聲成分的二階嵌套陣列信號模型，提出了改進的典范分解的嵌套陣列DOA估計算法，并對該算法進行了仿真，得到了相同信噪比和快拍數情況下，該算法和直接MUSIC、空間平滑算法，以及克拉美羅下界[14]的均方根誤差對比圖。仿真表明改進的典范分解算法適用于含噪條件下的二階嵌套陣列DOA估計，且與參考文獻[5，8]中提到的算法相比，具有更好的DOA估計性能。

2 二階嵌套陣列信號模型

N陣元的二階嵌套陣列結構如圖1所示，由兩個均勻線陣連接而成：第一階被稱為內部均勻線陣（inner ULA），含有N1個陣元且陣元間隔為d=λ/ 2 ，λ為接收信號波長；第二階被稱為外部均勻線陣（outer ULA）含有N2個陣元且陣元間隔為(N1+1)d。當N為偶數時，N1=N2=N/2；當N為奇數時N1=(N- 1 )/2，N2=(N+ 1 )/2。假設有R個窄帶互不相關的信號從方向為θi, (i=1 ,2,… ,R)的遠場空間入射該二階嵌套陣列，其K次接收快拍矩陣X∈CN×K表示為：

圖1 二階嵌套陣列結構

其中，矩陣A∈CN×R為入射信號的導向量矩陣A=[a(θ1),a(θ2) , … ,a(θR)]，向量a(θi)表示入射角度為θi的信號的導向量，向量si表示其K次采樣，矩陣S=[s1,s2,… ,sR]T∈CR×K，上標符號T表示矩陣轉秩，矩陣Wn∈CN×K表示所有陣元K次快拍接收的噪聲分量，噪聲為高斯白噪聲滿足正態分布，均值為0，方差為。

3 基于改進典范分解的嵌套陣列 DOA估計算法

3.1 典范分解算法基本原理

根據參考文獻[13]，假設K次接收快拍不含噪聲成分，對其做奇異值分解（singular value decomposition，SVD）得到上標符號H表示矩陣共軛轉秩，是有奇異值從大到小排列而成的對角矩陣，分別為與各個奇異值對應的左特征向量矩陣和右特征向量矩陣。一定存在 矩 陣使 得對 矩 陣做雙線性譜圖函數映射得到矩陣，雙線性譜圖函數映射目的是通過線性映射改變矩陣元素使其滿足典范分解的約束條件[12]，求解矩陣的核得到矩陣即求解矩陣廣義特征值，即上標符號-1表示矩陣求逆，故導向量矩陣A=，利用求解得到的導向量矩陣估計入射信號角度。

典范分解算法直接對接收快拍做SVD和雙線性譜圖函數映射，實質上是對接收數據做一階矩處理，在接收快拍不含噪聲成分時，矩陣Q不含有噪聲成分，為非滿秩矩陣，它的核存在非零解，但是當接收快拍中含有噪聲成分時，即接收快拍為X=AS+Wn而不再是=AS，一階距處理無法將信號和噪聲分離，矩陣含有噪聲成分，利用直接求核的方法得到的矩陣只含0元素，無法估計導向量矩陣和入射角度。

3.2 改進典范分解與嵌套陣列

得到二階嵌套陣列的K次接收快拍X后對其做SVD得到：

其中，∑∈CK×K是有奇異值從大到小排列而成的對角矩陣，矩陣U∈CN×R，矩陣V∈CK×K分別為與各個奇異值對應的左特征向量矩陣和右特征向量矩陣。UR=U(:,1:R)，Σs∈CR×R為信號空間對應的奇異值對角矩陣，且ΣR=Σ( 1:R,1:R)。參考文獻[11，14]，一定存在矩陣F∈CR×R，使得URΣR=AFT。定義矩陣Y∈CN×R且Y=URΣR。

Λ(θ)=表示R×R的對角矩陣，被稱為階間旋轉因子矩陣，其維度由入射信號個數決定，對角元素由入射信號角度和階間陣元距離差共同決定，根據二階嵌套陣列結構特性可得Y(2)=Y(1)×Λ(θ)。定義矩陣表示去掉矩陣第一行得到的子矩陣，=Y(1)(2:N1,:) , 矩陣表示去掉矩陣最后一行得到的子矩陣，同理可得矩陣表示向量J的第n個元素J=[N1- 1N2- 1 ]，定義矩陣

其中，i,j∈ [1,… ,Jn]，k,l∈ [1 ,2]。

如第 3.1節所述，矩陣Q含有噪聲成分，如果按照典范分解的方法采用直接求核的方法得到的矩陣只含0元素，無法估計導向量矩陣和入射角度。參考文獻[16]利用信號和噪聲在二階統計量空間正交的特性，采用奇異值分解將矩陣Q的協防差矩陣分解為正交的信號子空間和噪聲子空間，再利用噪聲子空間得到矩陣Mr，過程如下：對QHQ做奇異值分解得到中R個最小的奇異值構成遞增排列集合σ=[σ1σ2...σR],即σ1≤σ2≤ … ≤σR，該集合中每個元素作為奇異值時所對應的右特征向量被定義為,將m(r)的元素下角標定義為

重新排列各元素得到矩陣：

得到矩陣Mr后，將求解矩陣F∈CR×R的問題轉化問聯合對角化的問題，即求解矩陣G-1=F：

當R=2時，即入射信號個數為兩個，聯合對角化問題退化為廣義特征值分解（generalized eigenvalue decomposition，GEVD）問題，當R≥3時，需用迭代算法求解數值解，可參考文獻[16-18]。求得矩陣F后，利用URΣR=AFT，求解得到導向量矩陣A：

定義導向量矩陣A的第i列列向量為ci A，a(θ)(1)=i可得，定義以下函數則入射角iθ的計算式為：

算法步驟歸納如下：

步驟 1輸入接收快拍X=AS+Wn,做奇異值分解得到Y=URΣR；

步驟 2對Y進行分組重排得到Yred(n,r),利用函數φr,s(n),得到矩陣QHQ；

步驟3對QHQ進行奇異值分解，對奇異值和右奇異值向量進行重排得到矩陣Mr；

步驟4對Mr進行聯合對角化，求得矩陣F；

步驟5利用URΣR=AFT，求得導向量矩陣A；

步驟6利用式（13）求得入射角度θi。

改進的典范分解算法與傳統的典范算法相比，主要改進體現在步驟1和步驟3上。步驟1中，傳統算法的接收快拍不含噪聲，而實際應用中噪聲是普遍存在的，所以改進算法的接收快拍中加入噪聲。步驟2對接收快拍進行了分組和線性映射，傳統算法得到的不含噪聲，故在步驟 3中對其直接求核可以得到含有可行解的，而改進算法得到的Q中含有噪聲，故在步驟（3）中對其協方差矩陣進行SVD，利用正交特性得到含有角度信息的Mr，剩下步驟兩種算法相同。

假設接收快拍不含噪聲，在步驟3中對的協方差矩陣做SVD，得到的噪聲子空間奇異值均為0，然后求解其對應的右特征向量矩陣，問題等效為=0與直接求核結果一樣。因此在無噪條件下改進算法與傳統算法性能一致。

綜上所述，改進算法與傳統算法計算復雜度基本相當，但是傳統算法只能應用于無噪條件，而改進算法在無噪和有噪條件下均適用。

4 仿真實驗

仿真實驗采用8陣元二階嵌套陣列，即內部均勻線陣含有N1=4個陣元，各陣元間隔為d，外部均勻線陣含有N2=4，各陣元間隔為5d，陣元分布集合為[d, 2d, 3d, 4d, 5d, 10d, 15d, 20d]。仿真在不同信噪比和快拍次數的情況下，分別進行200次蒙特卡洛實驗，對比MUSIC、空間平滑和本文所提的改進典范分解算法 DOA估計結果的均方根誤差，并以克拉美羅下界值作為參考，定量分析算法估計性能。

實驗中信噪比（SNR）、均方根誤差（ERMS）和克拉美羅下界值（ C RBL）計算式分別為:

其中，表示二階范數，R表示信號個數，MC表示蒙特卡洛實驗次數，本文中MC=2 00，表示第mc次實驗所得估計角集合。

4.1 仿真實驗1

3個功率相同的信號分別從30°、45°、60°入射8陣元二階的嵌套陣列，信源數R=3，入射集合快拍次數K=50，信噪比從0 dB變化到20 dB，3種算法的均方根誤差如圖2所示。

圖2中，8陣元二階嵌套陣列比8陣元均勻線陣列擁有更高的自由度，因此在相同信噪比和快拍次數的情況下，8陣元二階嵌套陣列的克拉美羅界明顯低于8陣元均勻線陣列，體現了相同條件下，二階嵌套陣列的性能優越性。對于直接MUSIC、空間平滑和改進典范分解算法的均方根誤差，信噪比小于2 dB 時，3種算法均方根誤差接近；而信噪比大于2 dB時，本文所提的改進典范分解算法的均方根誤差明顯低于其余兩種算法；當信噪比大于10 dB時，改進典范分解算法對8陣元二階嵌套陣列的均方根誤差比直接 MUSIC算法小一個數量級，且低于8陣元均勻線陣的克拉美羅界。圖2驗證了改進典范分解算法對嵌套陣列的可行性和體現出性能優越性。

圖2 均方根誤差隨信噪比變化曲線

4.2 仿真實驗2

采用仿真1中的入射角度分別為30°、45°、60°的3個信號信噪比為10 dB，快拍次數從5變化到500，均方根誤差如圖3所示。

圖3 均方根誤差隨快拍數變化曲線

圖3中，當快拍次數小于50時，3種算法的均方根誤差受快拍次數變化影響較大，當快拍次數大于50時，變化趨勢趨于平緩，而實際情況中快拍數常常高于50，在這一區間內3種算法的均方根誤差受快拍數影響不大。仿真結果表明當信噪比為10 dB時，在相同快拍次數的情況下，改進典范分解算法均方根誤差小于直接 MUSIC算法和空間平滑算法，表明所提算法DOA估計精度高于另外兩種算法。

4.3 仿真實驗3

該仿真參考了文獻[19]所采用的在相同信噪比和快拍數的前提下，通過比較運行時間來對比運算復雜度的方式，使用同一計算機和同一軟件（實驗環境：CPU 2.5 GHz，8 GB RAM，MATLAB R2014a，Windows 8 x64）對本文所提算法、直接MUSIC算法和空間平滑算法，在快拍數為200和信噪比為10 dB（參考圖2和圖3，角度域搜索步長取0.05°）的情況下，分別做100次、200次和500次蒙特卡洛分析得到3種算法的運行時間，結果見表1。從表1可以看出，本文所提算法不需要峰值搜索，復雜度明顯低于直接MUSIC算法和空間平滑算法。

表1 算法運行時間比較

5 結束語

本文結合二階嵌套陣列含噪信號模型，首次提出了基于改進典范分解的嵌套陣列 DOA估計算法。該算法無需進行峰值搜索，利用局部平移不變特性，通過張量構造和典范分解求解得到導向量矩陣并求得入射角度。該算法利用SVD解將矩陣Q的協防差矩陣分解為正交的信號子空間和噪聲子空間,再利用噪聲子空間得到矩陣r M，解決了含噪條件下典范分解算法對嵌套陣列 DOA估計的問題，且在相同信噪比、快拍數情況下比直接MUSIC算法、空間平滑算法有更好的估計性能和更少的運算時間。