張 昊，吳晨曦

(1.安徽城市管理職業(yè)學(xué)院，安徽 合肥 230601；2.安徽省制約重點(diǎn)實(shí)驗(yàn)室，安徽 合肥 230037)

0 引言

波達(dá)方向估計(jì)(Direction of Arrival Estimation，DOA)是陣列信號(hào)處理領(lǐng)域研究的重要內(nèi)容之一，在雷達(dá)、通信、電子對(duì)抗等諸多領(lǐng)域都有著廣泛的應(yīng)用。近幾十年來，一系列高分辨方法相繼被提出，但現(xiàn)有大部分方法僅適用于超定條件(陣元數(shù)多于信號(hào)數(shù))， 隨著現(xiàn)代電磁環(huán)境的日益復(fù)雜，欠定DOA估計(jì)問題越來越常見，如何利用較少的陣元數(shù)對(duì)更多的信號(hào)進(jìn)行估計(jì)是目前研究的難點(diǎn)問題。

互質(zhì)陣列[1]的提出為解決欠定DOA估計(jì)問題提供了一個(gè)新的思路。由于其具有易于構(gòu)造、陣列擴(kuò)展性好、物理陣元和虛擬陣元具有解析表達(dá)式等優(yōu)點(diǎn)受到國內(nèi)外學(xué)者的廣泛關(guān)注[2-8]?，F(xiàn)有成果可以分為兩大類：一是基于子空間類方法[2-4]，該類方法的不足在于只能夠利用差聯(lián)合陣列中連續(xù)虛擬陣元部分而舍棄了非連續(xù)部分的虛擬陣元，導(dǎo)致陣列的虛擬孔徑存在一定的損失同時(shí)估計(jì)精度受到掃描網(wǎng)格的影響；二是基于稀疏重構(gòu)的方法[5-8]，該類方法雖然克服子空間方法的不足，但其前提條件是所有目標(biāo)準(zhǔn)確地位于預(yù)定的字典網(wǎng)格上，實(shí)際環(huán)境中并不一定滿足即存在模型失配問題，此時(shí)估計(jì)性能會(huì)嚴(yán)重下降。雖然字典網(wǎng)格越密集，DOA信息與字典網(wǎng)格之間的誤差越小，但計(jì)算量會(huì)隨之增加同時(shí)過于密集的網(wǎng)格會(huì)造成字典原子之間的相關(guān)性增強(qiáng)，影響重構(gòu)的性能。

陣列自由度與其相對(duì)應(yīng)的差聯(lián)合陣列中的虛擬陣元有關(guān)，而互質(zhì)陣列的差聯(lián)合陣列中由于存在空洞部分，若能夠?qū)⒖臻g部分進(jìn)行有效填充，進(jìn)而得到一個(gè)虛擬陣元更多、孔徑更大的均勻陣列，從而提高可估計(jì)信號(hào)數(shù)和估計(jì)精度[9-11]。文獻(xiàn)[9]通過陣列運(yùn)動(dòng)的方式來獲得與空缺項(xiàng)對(duì)應(yīng)的測(cè)量值，但該方法對(duì)運(yùn)動(dòng)狀態(tài)有嚴(yán)格的要求，在實(shí)際應(yīng)用中受限。文獻(xiàn)[10]通過引入額外的工作頻率來實(shí)現(xiàn)對(duì)空洞的填充，然而額外頻率的使用會(huì)增大系統(tǒng)的工作帶寬且會(huì)引起信號(hào)源特性起伏，在一定程度上增加了系統(tǒng)復(fù)雜性。文獻(xiàn)[11]提出一種基于稀疏重構(gòu)的外推技術(shù)，但該方法的估計(jì)性能同樣受到模型失配的影響[12]。針對(duì)上述問題，本文提出基于協(xié)方差矩陣重構(gòu)的互質(zhì)陣列DOA估計(jì)方法。

1 信號(hào)模型

互質(zhì)陣列由兩個(gè)不同的均勻子陣組成，如圖1所示，其中子陣1包含有N個(gè)物理陣元，分別位于{Mnd,0≤n≤N-1}，子陣2包含有2M-1個(gè)物理陣元，分別位于{Nmd,1≤m≤2M-1}，M,N是兩個(gè)互質(zhì)的整數(shù)，d=λ/2，λ為入射信號(hào)波長。

圖1 互質(zhì)陣列結(jié)構(gòu)圖Fig.1 Structure of coprime array

假設(shè)K個(gè)遠(yuǎn)場(chǎng)窄帶平面波信號(hào)分別以[θ1,θ2,…,θK]入射到互質(zhì)陣列上，t時(shí)刻陣列接收數(shù)據(jù)可表示為：

X(t)=AS(t)+n(t)

(1)

式(1)中，X(t)=[x1(t),x2(t),…,x2M+N-1(t)]T為陣列接收數(shù)據(jù)向量，S(t)=[s1(t),…,sK(t)]T為空間信號(hào)向量，n(t)=[n1(t),n2(t),…,n2M+N-1(t)]T為噪聲向量，A=[a(θ1),…,a(θK)]為(2M+N-1)×K維陣列流型矩陣，其中a(θi)為入射角θi的導(dǎo)向矢量。

假設(shè)各入射信號(hào)之間互不相關(guān)，且噪聲為高斯白噪聲且與信號(hào)相互獨(dú)立。則陣列接收數(shù)據(jù)的協(xié)方差矩陣可表示為：

RX=E[X(t)XH(t)]=ARsAH+n

(2)

在實(shí)際應(yīng)用中，由于快拍數(shù)L有限，RX由采樣協(xié)方差矩陣進(jìn)行等效，即

(3)

2 基于協(xié)方差矩陣重構(gòu)的DOA估計(jì)方法

首先，給出差聯(lián)合陣列(difference coarray)的定義，定義如下集合

D={xi-xj},xi,xj∈C

(4)

式(4)中，xi為第i個(gè)陣元位置，C為所有陣元位置集合。D表示陣列中所有陣元位置差值構(gòu)成的集合。由于集合D中存在冗余元素，定義集合D中所有不相同的元素組成的集合Du為差聯(lián)合陣列。

圖1所示的互質(zhì)陣列的差聯(lián)合陣列如圖2所示，虛擬陣元的位置可具體表示為{±(Mnd-Nmd),0≤n≤N-1,0≤m≤2M-1}。

圖2 差聯(lián)合陣列示意圖Fig.2 Structure of difference coarray

根據(jù)虛擬陣元分布可知，差聯(lián)合陣列可以被劃分為一個(gè)均勻陣列和兩個(gè)非均勻陣列。其中，均勻陣列包含有2MN+2M-1個(gè)虛擬陣元以陣元間距d均勻分布且關(guān)于原點(diǎn)對(duì)稱。而非均勻陣列由于存在空洞從而造成虛擬陣元呈現(xiàn)非均勻分布，兩個(gè)非均勻陣列同樣關(guān)于原點(diǎn)對(duì)稱。進(jìn)一步可知，非均勻陣列可再細(xì)分為單個(gè)空洞部分和多個(gè)空洞部分?，F(xiàn)有的大部分方法僅利用了連續(xù)部分的虛擬陣元，而未得到有效利用的虛擬陣元數(shù)為MN+M+N-1，因此，若能將差聯(lián)合陣列存在的空洞進(jìn)行有效填充，則可顯著提高可利用的虛擬陣元數(shù)，從而獲得更大的陣列自由度，實(shí)現(xiàn)對(duì)更多信號(hào)的DOA估計(jì)。

實(shí)際上差聯(lián)合陣列與陣列所能得到的波程差是一一對(duì)應(yīng)的。對(duì)于圖1所示的互質(zhì)陣列，陣列接收數(shù)據(jù)協(xié)方差矩陣RX具體可表示為：

(5)

在實(shí)際應(yīng)用中，由于受到信噪比和快拍數(shù)有限的影響，式協(xié)方差矩陣中同一波程差對(duì)應(yīng)的元素并不完全相等，因此，為了充分利用樣本信息提高估計(jì)的精度，對(duì)相同波程差的元素進(jìn)行求平均運(yùn)算作為該波程差對(duì)應(yīng)的輸出值，即

(6)

式(6)中，Ri(du)表示對(duì)應(yīng)波程差R(du)的第i個(gè)元素。

進(jìn)一步根據(jù)存在的波程差元素對(duì)協(xié)方差矩陣RX進(jìn)行擴(kuò)展，構(gòu)成一個(gè)(2MN-N+1)×(2MN-N+1)維的Toeplitz矩陣RT：

(7)

矩陣RT可看作一個(gè)數(shù)據(jù)缺失的協(xié)方差矩陣，同時(shí)由于入射信號(hào)的稀疏性，因此，RT是低秩的。使得利用低秩矩陣重構(gòu)理論[13]對(duì)矩陣中的零元素進(jìn)行填充成為可能，具體可通過求解如下的優(yōu)化問題：

(8)

式(8)中，R表示待重構(gòu)的目標(biāo)矩陣，η表示誤差門限，該值的選取對(duì)算法性能的具有重要的影響。P是映射矩陣，一旦互質(zhì)陣列結(jié)構(gòu)確定，投影映射矩陣隨之確定，具體可表示為：

(9)

然而，矩陣秩函數(shù)是非連續(xù)、非凸的，直接求解矩陣秩最小化問題是一個(gè)NP-hard問題。這里，利用跡范數(shù)最小化來對(duì)秩范數(shù)最小化進(jìn)行凸松弛表示。因此，式(8)進(jìn)一步可表示為：

(10)

vec(E)～AsN(0,W)

(11)

(12)

進(jìn)一步可表示為：

(13)

式(13)中，Asχ2((2MN-N)2)表示自由度為(2MN-N)2的近似卡方分布。進(jìn)一步在式(13)引入?yún)?shù)μ：

(14)

式(14)中，μ是一個(gè)決定噪聲容錯(cuò)的參數(shù)，(在Matlab軟件中可利用函數(shù)chi2inv(1-ρ,(2MN-N)2)計(jì)算)，μ與置信關(guān)聯(lián)ρ相關(guān)聯(lián)。為了獲得較高的概率1-ρ，ρ的取值一般很小。例如，ρ=0.001。

通過將式(14)代入式(10)，DOA估計(jì)式可表示為：

mintr(R)

(15)

式(15)是一個(gè)凸優(yōu)化問題，可通過CVX[15]優(yōu)化工具箱進(jìn)行求解，得到最優(yōu)估計(jì)值R。

進(jìn)一步對(duì)R進(jìn)行特征分解，即

(16)

式(6)中，Us是大特征值對(duì)應(yīng)的特征向量所構(gòu)成的(2MN-N+1)×K維信號(hào)空間，UN是小特征值對(duì)應(yīng)的特征向量所構(gòu)成的(2MN-N+1)×(2MN-N+1-K)維噪聲空間，Σs是大特征值組成的K×K維對(duì)角陣，ΣN為小特征值組成的(2MN-N+1-K)×(2MN-N+1-K)維對(duì)角陣。

為了避免角度搜索的復(fù)雜性，令P(z)=[1,z,…,z2MN-N]T，則DOA的求根多項(xiàng)式為：

(17)

(18)

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

本節(jié)中將通過仿真實(shí)驗(yàn)對(duì)所提方法的性能進(jìn)行驗(yàn)證，并與SS-MUSIC方法、Lasso方法以及文獻(xiàn)[11]中的EX-Lasso方法進(jìn)行比較，其中，SS-MUSIC方法的角度搜索間隔為0.1°，Lasso方法和EX-Lasso方法的完備字典間隔為0.1°。以互質(zhì)陣列作為接收陣列，最小陣元間距d為半波長。

實(shí)驗(yàn)1 為驗(yàn)證本文方法在欠定條件下的估計(jì)性能，假設(shè)9個(gè)等功率遠(yuǎn)場(chǎng)窄帶信號(hào)入射到互質(zhì)陣列上，陣元數(shù)為6，其中，M=2，N=3，入射角度分別為[-70.12°，-50.34°，-30.02°，-10.25°，10.65°，30.45°，50.12°，60.34°，80.23°]，信噪比SNR=10 dB，快拍數(shù)L=1 024，實(shí)驗(yàn)結(jié)果如圖3所示。

圖3 四種方法可行性比較Fig.3 Feasibility comparison of four methods

由圖3的實(shí)驗(yàn)結(jié)果可知，本文方法利用6個(gè)陣元準(zhǔn)確實(shí)現(xiàn)了對(duì)9個(gè)入射信號(hào)DOA的估計(jì)，而MUSIC方法理論上最多可估計(jì)信號(hào)數(shù)為7，Lasso方法的可估計(jì)信號(hào)數(shù)為8均小于入射信號(hào)數(shù)，均無法實(shí)現(xiàn)對(duì)所有信號(hào)DOA的準(zhǔn)確估計(jì)。EX-Lasso方法雖然理論上能夠?qū)崿F(xiàn)對(duì)全部信號(hào)的DOA估計(jì)，但在重構(gòu)過程中會(huì)受到模型失配的影響，從而對(duì)部分入射信號(hào)的DOA估計(jì)效果不佳。

實(shí)驗(yàn)2 為驗(yàn)證本文方法的估計(jì)精度與信噪比、快拍數(shù)之間的關(guān)系。假設(shè)2個(gè)等功率遠(yuǎn)場(chǎng)窄帶信號(hào)入射到互質(zhì)陣列上，陣元數(shù)為6，其中M=2，N=3，入射角度分別為[30.13°，60.55°]，分析信噪比、快拍數(shù)對(duì)四種方法估計(jì)精度的影響，每個(gè)信噪比、快拍數(shù)下進(jìn)行500次蒙特卡羅實(shí)驗(yàn)。圖4為快拍數(shù)L=1 024，信噪比SNR從-5～15 dB，步長為5 dB時(shí)，角度均方根誤差隨信噪比變化曲線；圖5為信噪比SNR=10 dB，快拍數(shù)L從500～900，步長為100時(shí)，角度均方根誤差隨快拍數(shù)變化曲線。

圖4 角度均方根誤差隨信噪比變化曲線Fig.4 RMSE versus SNR

圖5 角度均方根誤差隨快拍數(shù)變化曲線Fig.5 RMSE versus snapshot number

由圖4和圖5的實(shí)驗(yàn)結(jié)果可知，在一定的信噪比和快拍數(shù)范圍內(nèi)，四種方法的估計(jì)均方根誤差均隨著信噪比的增加和快拍數(shù)的增大而減少。在相同的信噪比和快拍數(shù)條件下，本文方根的估計(jì)精度要明顯高于其他三種方法，尤其是在低信噪比和少快拍數(shù)時(shí)，從而驗(yàn)證了本文方法具有較好的估計(jì)性能。這是由于一方面，本文方法對(duì)互質(zhì)陣列的差聯(lián)合陣列空洞部分的有效填充，提高了可利用的陣列自由度；另一方面，SS-MUSIC方法的估計(jì)精度受到角度搜索間隔的影響，Lasso方法和EX-Lasso方法在重構(gòu)過程中均受到模型失配的影響，而本文文法有效避免了該問題，從而具有更高的估計(jì)精度。

實(shí)驗(yàn)3 為驗(yàn)證本文方法對(duì)相近信號(hào)的分辨能力，假設(shè)2個(gè)相近的等功率遠(yuǎn)場(chǎng)窄帶信號(hào)入射到互質(zhì)陣列上，陣元數(shù)為6，其中M=2，N=3，入射角度分別為[30.13°、32.26°]，信噪比SNR=5 dB，快拍數(shù)L=500。將四種方法的分辨性能進(jìn)行比較，圖6為四種方法具體的分辨結(jié)果。

圖6 四種方法的分辨性能比較Fig.6 Comparison of the resolution of four methods

由圖6的實(shí)驗(yàn)結(jié)果可知，當(dāng)2個(gè)入射信號(hào)較近時(shí)，本文方法仍然能夠準(zhǔn)確估計(jì)出入射信號(hào)的DOA信息，而其他三種方法此時(shí)失效，體現(xiàn)了本文方法對(duì)相近信號(hào)較好的分辨能力。

4 結(jié)論

本文提出了基于協(xié)方差矩陣重構(gòu)的互質(zhì)陣列DOA估計(jì)方法，該方法首先根據(jù)差聯(lián)合陣列與波程差一一對(duì)應(yīng)的關(guān)系，將陣列協(xié)方差矩陣擴(kuò)展為一個(gè)Toeplitz矩陣，進(jìn)一步通過求解跡范數(shù)最小化問題實(shí)現(xiàn)對(duì)矩陣中零元素的恢復(fù)，最后利用求根MUSIC方法實(shí)現(xiàn)對(duì)DOA的估計(jì)。仿真實(shí)驗(yàn)結(jié)果表明，該方法對(duì)互質(zhì)陣列的差聯(lián)合陣列中空洞部分能夠有效進(jìn)行填充，從而獲得更高的陣列自由度。在提高可估計(jì)信號(hào)數(shù)的同時(shí)，有效避免了傳統(tǒng)稀疏重構(gòu)中模型失配對(duì)估計(jì)性能的影響，提高了估計(jì)精度和分辨性能。