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

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

1 引言

在陣列信號(hào)處理中，更高的自由度（degree of freedoms，DOF）往往能夠帶來(lái)更高的輸出信噪比、更低的旁瓣、更低的克拉美羅下界（Cramer-Rao lower bound，CRLB）等。與均勻線陣（uniformlinear array，ULA）相比，非均勻線陣（nonuniform linear array，NLA）因?yàn)槠湎∈杼匦阅軌蛟谑褂孟嗤囋獢?shù)的情況下，獲得更高的 DOF[1]。常見(jiàn)的非均勻線陣有最小不冗余陣（minimum redundancy array，MRA）[2]、互質(zhì)陣列（co-prime array）[3-4]以及嵌套陣列（nested array）[5]。

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

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

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

2 二階嵌套陣列信號(hào)模型

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

圖1 二階嵌套陣列結(jié)構(gòu)

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

3 基于改進(jìn)典范分解的嵌套陣列 DOA估計(jì)算法

3.1 典范分解算法基本原理

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

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

3.2 改進(jìn)典范分解與嵌套陣列

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

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

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

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

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

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

得到矩陣Mr后，將求解矩陣F∈CR×R的問(wèn)題轉(zhuǎn)化問(wèn)聯(lián)合對(duì)角化的問(wèn)題，即求解矩陣G-1=F：

當(dāng)R=2時(shí)，即入射信號(hào)個(gè)數(shù)為兩個(gè)，聯(lián)合對(duì)角化問(wèn)題退化為廣義特征值分解（generalized eigenvalue decomposition，GEVD）問(wèn)題，當(dāng)R≥3時(shí)，需用迭代算法求解數(shù)值解，可參考文獻(xiàn)[16-18]。求得矩陣F后，利用URΣR=AFT，求解得到導(dǎo)向量矩陣A：

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

算法步驟歸納如下：

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

步驟 2對(duì)Y進(jìn)行分組重排得到Y(jié)red(n,r),利用函數(shù)φr,s(n),得到矩陣QHQ；

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

步驟4對(duì)Mr進(jìn)行聯(lián)合對(duì)角化，求得矩陣F；

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

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

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

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

綜上所述，改進(jìn)算法與傳統(tǒng)算法計(jì)算復(fù)雜度基本相當(dāng)，但是傳統(tǒng)算法只能應(yīng)用于無(wú)噪條件，而改進(jìn)算法在無(wú)噪和有噪條件下均適用。

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

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

實(shí)驗(yàn)中信噪比（SNR）、均方根誤差（ERMS）和克拉美羅下界值（ C RBL）計(jì)算式分別為:

其中，表示二階范數(shù)，R表示信號(hào)個(gè)數(shù)，MC表示蒙特卡洛實(shí)驗(yàn)次數(shù)，本文中MC=2 00，表示第mc次實(shí)驗(yàn)所得估計(jì)角集合。

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

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

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

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

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

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

圖3 均方根誤差隨快拍數(shù)變化曲線

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

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

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

表1 算法運(yùn)行時(shí)間比較

5 結(jié)束語(yǔ)

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