基于協方差矩陣重構的互質陣列DOA估計方法

張 昊，吳晨曦

(1.安徽城市管理職業學院，安徽 合肥 230601；2.安徽省制約重點實驗室，安徽 合肥 230037)

0 引言

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

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

陣列自由度與其相對應的差聯合陣列中的虛擬陣元有關，而互質陣列的差聯合陣列中由于存在空洞部分，若能夠將空間部分進行有效填充，進而得到一個虛擬陣元更多、孔徑更大的均勻陣列，從而提高可估計信號數和估計精度[9-11]。文獻[9]通過陣列運動的方式來獲得與空缺項對應的測量值，但該方法對運動狀態有嚴格的要求，在實際應用中受限。文獻[10]通過引入額外的工作頻率來實現對空洞的填充，然而額外頻率的使用會增大系統的工作帶寬且會引起信號源特性起伏，在一定程度上增加了系統復雜性。文獻[11]提出一種基于稀疏重構的外推技術，但該方法的估計性能同樣受到模型失配的影響[12]。針對上述問題，本文提出基于協方差矩陣重構的互質陣列DOA估計方法。

1 信號模型

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

圖1 互質陣列結構圖Fig.1 Structure of coprime array

假設K個遠場窄帶平面波信號分別以[θ1,θ2,…,θK]入射到互質陣列上，t時刻陣列接收數據可表示為：

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

(1)

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

假設各入射信號之間互不相關，且噪聲為高斯白噪聲且與信號相互獨立。則陣列接收數據的協方差矩陣可表示為：

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

(2)

在實際應用中，由于快拍數L有限，RX由采樣協方差矩陣進行等效，即

(3)

2 基于協方差矩陣重構的DOA估計方法

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

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

(4)

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

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

圖2 差聯合陣列示意圖Fig.2 Structure of difference coarray

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

實際上差聯合陣列與陣列所能得到的波程差是一一對應的。對于圖1所示的互質陣列，陣列接收數據協方差矩陣RX具體可表示為：

(5)

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

(6)

式(6)中，Ri(du)表示對應波程差R(du)的第i個元素。

進一步根據存在的波程差元素對協方差矩陣RX進行擴展，構成一個(2MN-N+1)×(2MN-N+1)維的Toeplitz矩陣RT：

(7)

矩陣RT可看作一個數據缺失的協方差矩陣，同時由于入射信號的稀疏性，因此，RT是低秩的。使得利用低秩矩陣重構理論[13]對矩陣中的零元素進行填充成為可能，具體可通過求解如下的優化問題：

(8)

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

(9)

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

(10)

vec(E)～AsN(0,W)

(11)

(12)

進一步可表示為：

(13)

式(13)中，Asχ2((2MN-N)2)表示自由度為(2MN-N)2的近似卡方分布。進一步在式(13)引入參數μ：

(14)

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

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

mintr(R)

(15)

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

進一步對R進行特征分解，即

(16)

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

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

(17)

(18)

3 仿真實驗

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

實驗1 為驗證本文方法在欠定條件下的估計性能，假設9個等功率遠場窄帶信號入射到互質陣列上，陣元數為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，快拍數L=1 024，實驗結果如圖3所示。

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

由圖3的實驗結果可知，本文方法利用6個陣元準確實現了對9個入射信號DOA的估計，而MUSIC方法理論上最多可估計信號數為7，Lasso方法的可估計信號數為8均小于入射信號數，均無法實現對所有信號DOA的準確估計。EX-Lasso方法雖然理論上能夠實現對全部信號的DOA估計，但在重構過程中會受到模型失配的影響，從而對部分入射信號的DOA估計效果不佳。

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

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

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

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

實驗3 為驗證本文方法對相近信號的分辨能力，假設2個相近的等功率遠場窄帶信號入射到互質陣列上，陣元數為6，其中M=2，N=3，入射角度分別為[30.13°、32.26°]，信噪比SNR=5 dB，快拍數L=500。將四種方法的分辨性能進行比較，圖6為四種方法具體的分辨結果。

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

由圖6的實驗結果可知，當2個入射信號較近時，本文方法仍然能夠準確估計出入射信號的DOA信息，而其他三種方法此時失效，體現了本文方法對相近信號較好的分辨能力。

4 結論

本文提出了基于協方差矩陣重構的互質陣列DOA估計方法，該方法首先根據差聯合陣列與波程差一一對應的關系，將陣列協方差矩陣擴展為一個Toeplitz矩陣，進一步通過求解跡范數最小化問題實現對矩陣中零元素的恢復，最后利用求根MUSIC方法實現對DOA的估計。仿真實驗結果表明，該方法對互質陣列的差聯合陣列中空洞部分能夠有效進行填充，從而獲得更高的陣列自由度。在提高可估計信號數的同時，有效避免了傳統稀疏重構中模型失配對估計性能的影響，提高了估計精度和分辨性能。