基于協(xié)方差矩陣重構(gòu)的互質(zhì)陣列DOA估計(jì)

盤敏容， 蔣留兵2,3， 車 俐2,3， 姜 興

(1.桂林電子科技大學(xué)信息與通信學(xué)院， 廣西桂林 541004； 2.桂林電子科技大學(xué)廣西無(wú)線寬帶通信與信號(hào)處理重點(diǎn)實(shí)驗(yàn)室， 廣西桂林 541004； 3.桂林電子科技大學(xué)計(jì)算機(jī)與信息安全學(xué)院， 廣西桂林 541004)

0 引言

在陣列信號(hào)處理中，波達(dá)方向(Direction of Arrival,DOA)估計(jì)算法被廣泛應(yīng)用于雷達(dá)、聲吶和無(wú)線通信等領(lǐng)域，基于均勻線陣(Uniform Linear Array，ULA)的DOA估計(jì)技術(shù)成熟且陣列結(jié)構(gòu)簡(jiǎn)單，是最常用的陣列DOA估計(jì)方法[1-2]。近年來(lái)，由于互質(zhì)陣列在相同陣元數(shù)下較均勻陣列有更大的孔徑，由于其出色的性能引起了廣泛的關(guān)注[3]。然而，由于互質(zhì)陣列的非均勻性，該陣列相應(yīng)的DOA估計(jì)算法較均勻線陣更難實(shí)現(xiàn)，故互質(zhì)陣列DOA估計(jì)算法具有一定研究意義。將互質(zhì)陣列通過(guò)某種方式拓展成虛擬陣列，虛擬陣列陣元數(shù)比實(shí)際的陣元數(shù)多，在一定程度上可以增加自由度(Degrees-of-Freedom，DOF)。由互質(zhì)陣拓展得到的虛擬陣列由于缺失陣元而不是線性均勻陣列，運(yùn)用傳統(tǒng)的DOA估計(jì)方法無(wú)法得出準(zhǔn)確的波達(dá)方向。為了解決這個(gè)問(wèn)題，利用虛擬陣列的最大連續(xù)陣元的信息進(jìn)行空域平滑的DOA估計(jì)，但該方法會(huì)損失不連續(xù)陣元的信息，影響DOA估計(jì)的自由度和分辨率[4]。文獻(xiàn)[5]通過(guò)將組成互質(zhì)陣中的較小均勻線陣的個(gè)數(shù)加倍，用虛擬陣列的最大連續(xù)陣元的信息進(jìn)行空域平滑的DOA估計(jì)提高了自由度和分辨率，但實(shí)際陣元數(shù)的增加會(huì)占用更多的空間資源。本文提出的算法在保持實(shí)際陣元數(shù)不變的情況下，通過(guò)陣列內(nèi)插的方式將虛擬陣列缺失的陣元補(bǔ)齊從而達(dá)到線性均勻的虛擬陣列，重構(gòu)協(xié)方差矩陣并優(yōu)化相關(guān)的元素進(jìn)行DOA估計(jì)，提高了自由度和分辨率。

1 信號(hào)模型

圖1給出了互質(zhì)陣的陣列結(jié)構(gòu)。該陣列由圖1(a)中的兩個(gè)稀疏均勻線陣(ULA)嵌套構(gòu)成，記為子陣1和子陣2，互質(zhì)陣中子陣1和子陣2的第1個(gè)陣元相互重合。子陣1有M個(gè)陣元，陣元間距為Nd，子陣2有N個(gè)陣元，陣元間距為Md，M和N為互質(zhì)的整數(shù)，d=λ/2，λ一般為半波長(zhǎng)。由于M和N的互質(zhì)性，除了第一個(gè)陣元，該互質(zhì)陣列其他陣元不會(huì)相互疊加，陣元數(shù)為M+N-1。

A(θ)s(t)+n(t)

(1)

計(jì)算互質(zhì)陣列接收信號(hào)x(t)的協(xié)方差矩陣如式(2):

Rx=E(x(t)xH(t) )=

(2)

由于Rx在實(shí)際應(yīng)用中無(wú)法直接獲得，在經(jīng)典DOA估計(jì)中，其值通常由樣本協(xié)方差矩陣代替為[6]

(3)

式中,T表示快拍數(shù)。

2 互質(zhì)陣列DOA估計(jì)

目前利用互質(zhì)陣列進(jìn)行DOA估計(jì)的方法主要有兩種[7]，一是利用互質(zhì)陣的子陣列進(jìn)行DOA估計(jì)，聯(lián)合兩個(gè)子陣列的估計(jì)結(jié)果進(jìn)行角度的解模糊，但該方法對(duì)信源分辨的個(gè)數(shù)小于子陣列的個(gè)數(shù)；二是對(duì)互質(zhì)陣列接收信號(hào)的協(xié)方差矩陣矢量化后構(gòu)造虛擬陣元，然后利用該虛擬陣列進(jìn)行DOA估計(jì)。為了充分利用互質(zhì)陣列的信息，更大地提高其自由度，本文采用方法二進(jìn)行分析。

2.1 互質(zhì)陣列的虛擬陣元拓展方式

(4)

Ldiff1={±(Mn-Nm)|n=0,1,…,N-1,

m=0,1,…,M-1}

(5)

該差分集對(duì)應(yīng)的矢量

(6)

式中，C1為位于Ldiff1d的互質(zhì)虛擬陣列的導(dǎo)向矢量構(gòu)成的矩陣，i1表示iv中取出重復(fù)陣元后的陣列對(duì)應(yīng)的元素。

圖2 互質(zhì)陣列以及拓展的虛擬陣列

現(xiàn)以子陣1陣元數(shù)M=3,子陣2陣元數(shù)N=5為例，對(duì)互質(zhì)陣列的虛擬陣元拓展以及虛擬陣列的DOA估計(jì)進(jìn)行分析，則實(shí)際陣元的位置為{0,3,5,6,9,10，12}。由此集合得到具有重復(fù)元素的差分集L1為{0,-3,-5,-6,-9,-10,-12,3,0,-2,-3,…,1,0,-2,12,9,7,6,3,2,0}，重復(fù)陣元位置的元素取均值并進(jìn)行重排處理得到差分集L2為{-12,-10,-9,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,9,10,12}，如圖2(b)所示。雖然圖2(b)中的陣列相對(duì)于互質(zhì)陣列陣元數(shù)增加，但在位置{-11,-8,8,11}處存在虛擬陣元的缺失，導(dǎo)致該虛擬陣列是一個(gè)非均勻陣列。

2.2 最大連續(xù)虛擬陣列空域平滑DOA估計(jì)

如圖3所示，將該陣列分解成L=(|L3|+1)/2個(gè)相互重疊的均勻子線陣，|·|表示集合的基數(shù)，每一個(gè)均勻子線陣包含L個(gè)連續(xù)虛擬陣元，v2,p，l=1,2,…,L分別為L(zhǎng)個(gè)子陣的矢量信息。進(jìn)行空域平滑得到矩陣：

(7)

最后利用子空間類算法對(duì)矩陣Rv進(jìn)行特征分解，張成對(duì)應(yīng)的子空間，從而實(shí)現(xiàn)DOA估計(jì)。

圖3 最大連續(xù)虛擬陣元的空域平滑

2.3 基于虛擬陣列內(nèi)插的協(xié)方差矩陣重構(gòu)DOA估計(jì)

利用最大連續(xù)虛擬陣元的矢量信息進(jìn)行DOA估計(jì)的方法會(huì)造成部分接收信號(hào)信息的丟失，從而影響后續(xù)DOA估計(jì)的精度并降低自由度。為了充分利用虛擬陣元的信息，我們引入了陣列內(nèi)插的思想來(lái)構(gòu)造一個(gè)包含L1中所有虛擬陣元的均勻陣列，形成如圖4所示的虛擬陣列L4。建立L4的向量化協(xié)方差矩陣v3和L1向量化協(xié)方差矩陣v1之間的關(guān)系，原虛擬陣元的等效信號(hào)不變，插入陣元接收到的信號(hào)為0[9]:

(8)

同理，v3是一個(gè)秩為1的向量，無(wú)法直接分辨多個(gè)信號(hào)源，需要通過(guò)圖4空域平滑來(lái)解決。由于v3關(guān)于原點(diǎn)對(duì)稱，利用式(9)簡(jiǎn)化初始的協(xié)方差矩陣，該矩陣包含了線性虛擬陣元中的全部信息，故省略空間平滑過(guò)程，即

R2=T(v3+)

(9)

(10)

式中，ε為協(xié)方差矩陣擬合誤差，rank(·)為矩陣的秩。

圖4 內(nèi)插虛擬陣列空域平滑

為簡(jiǎn)化式(10)中的優(yōu)化問(wèn)題，定義與矩陣R2相同維度的投影矩陣P，矩陣R2與P中的元素均與圖4的內(nèi)插均勻陣列的矢量元素對(duì)應(yīng)，但矩陣P對(duì)應(yīng)內(nèi)插陣元的元素值為0，其余為1。引入矩陣的跡，將式(10)中的優(yōu)化問(wèn)題轉(zhuǎn)化為凸優(yōu)化問(wèn)題[10]：

s.t.T(z)0

(11)

式中,τ為正則化參數(shù)，‖·‖F(xiàn)表示計(jì)算Frobenius范數(shù)。在計(jì)算以向量z為第一列的厄米特托普利茲矩陣時(shí)，由于z是由復(fù)數(shù)元素組成的向量，該向量的第一個(gè)元素?zé)o法進(jìn)行共軛拓展，故先將向量z的第一個(gè)元素定義為實(shí)數(shù)，并且在完成向量z的整體優(yōu)化后再單獨(dú)優(yōu)化向量z的第一個(gè)元素。最后利用子空間類算法對(duì)優(yōu)化重構(gòu)后的矩陣T(z)進(jìn)行特征分解，張成對(duì)應(yīng)的子空間，進(jìn)行譜峰搜索，得到波達(dá)方向。

2.4 算法步驟總結(jié)

根據(jù)上述分析過(guò)程，將基于虛擬陣列內(nèi)插的協(xié)方差矩陣重構(gòu)DOA估計(jì)算法的主要步驟歸納如下：

步驟1： 根據(jù)互質(zhì)陣模型以及預(yù)設(shè)入射角，求解接收信號(hào)的協(xié)方差矩陣，如式(3)并對(duì)協(xié)方差矩陣進(jìn)行向量化處理，如式(4)，初步構(gòu)造虛擬陣列模型，如式(6)；

步驟2： 對(duì)不連續(xù)的虛擬陣列進(jìn)行內(nèi)插，在構(gòu)造的均勻虛擬陣列的基礎(chǔ)上初始化協(xié)方差矩陣，如式(8)、式(9)；

步驟3： 利用初始化的協(xié)方差矩陣進(jìn)行凸優(yōu)化，初步重構(gòu)協(xié)方差矩陣，如式(11)，然后優(yōu)化矢量z的第一個(gè)元素以完成矩陣T(z)的重構(gòu)；

步驟4： 對(duì)重構(gòu)后的矩陣T(z)進(jìn)行特征值分解，依據(jù)信號(hào)子空間和噪聲子空間的正交性，進(jìn)行譜峰搜索，最終得到DOA估計(jì)結(jié)果。

3 實(shí)驗(yàn)和結(jié)果分析

為了驗(yàn)證算法的性能，分別對(duì)最大連續(xù)虛擬陣列(Maximum Continuous Virtual Array，MCVA)DOA估計(jì)算法和協(xié)方差矩陣重構(gòu)(Covariance Matrix Reconstruction，CMR)DOA估計(jì)算法進(jìn)行仿真對(duì)比分析?；谏鲜鲫嚵心Ｐ秃退惴ㄟM(jìn)行了仿真分析,其中兩組互質(zhì)的稀疏均勻陣列的陣元數(shù)分別為M=3,N=5，由這兩組陣列組成的非均勻陣列的陣元數(shù)為M+N-1=7，分別位于{0,3d,5d,6d,9d,10d,12d}，d=λ/2。正則化參數(shù)τ的取值為2.2×10-5，步驟3對(duì)矢量z第一個(gè)元素最后的優(yōu)化值為6.0。

實(shí)驗(yàn)1 不同信源數(shù)時(shí)算法的檢測(cè)性能仿真實(shí)驗(yàn)

為了驗(yàn)證兩種算法的檢測(cè)能力，設(shè)置仿真條件：信噪比SNR=5 dB,快拍數(shù)K=1 000,來(lái)波分布在-50°～50°之間，對(duì)信源數(shù)為K=7和K=9作了仿真。圖5給出在多目標(biāo)情況下利用MCVA算法和CMR算法計(jì)算出的空間譜，豎線為所設(shè)置的信號(hào)源入射角。當(dāng)M=3,n=5時(shí)，是式(7)中MCVA算法空域平滑計(jì)算出的矩陣Rv的最大秩為8，故該算法最多能分辨7個(gè)不同方向的信號(hào)源。由圖5可知，MCVA算法和CMR算法均能分辨7個(gè)不同方向信號(hào)源，且CMR算法可分辨9個(gè)不同方向信號(hào)源，故CMR算法具有更大的自由度。

圖5 不同信源數(shù)K下MCVA算法和CMR算法的空間譜

實(shí)驗(yàn)2 算法的分辨率性能仿真實(shí)驗(yàn)

分辨率表征了算法對(duì)相鄰信號(hào)源的譜峰的分辨能力。設(shè)置仿真條件：信噪比SNR=5 dB,快拍數(shù)K=1 000,研究?jī)煞N算法在相鄰信號(hào)源方向分別為θ1=-3°,θ2=3°，θ1=-2°,θ2=2°，θ1=-1.5°,θ2=1.5°的情況下的分辨情況，實(shí)驗(yàn)結(jié)果如圖6所示，豎線表示信號(hào)源入射方向。由圖6可知，MCVA算法在θ1=-2°,θ2=2°時(shí)無(wú)法分辨兩個(gè)目標(biāo)，在θ1=-3°,θ2=3°可分辨兩個(gè)目標(biāo)，而CMR算法在θ1=-1.5°,θ2=1.5°和θ1=-2°,θ2=2°均能分辨兩個(gè)目標(biāo)，故CMR算法的分辨率較高。

圖6 算法在相鄰目標(biāo)θ1,θ2下的空間譜

實(shí)驗(yàn)3 算法的估計(jì)精度仿真實(shí)驗(yàn)

為了衡量算法的檢測(cè)精度，定義均方根誤差(Root Mean Square Error,RMSE)的數(shù)學(xué)表達(dá)式如式(12)：

(12)

(a) 不同信噪比

(b) 不同快拍數(shù)圖7 算法在信源數(shù)K=2下的均方根誤差

(a) 信源數(shù)K=7

(b) 信源數(shù)K=9圖8 多信源數(shù)下算法在不同快拍數(shù)下的性能

如圖8(a)所示為信源數(shù)K=7,方向?yàn)棣? [-50°,-35°,-15°,0°,10°,30°,50°]，信噪比SNR=5 dB，Monte-Carlo試驗(yàn)次數(shù)為P=30下MCVA算法和CMR算法在不同快拍數(shù)下的RMSE的實(shí)驗(yàn)結(jié)果；如圖8(b)所示為信源數(shù)K=9,方向?yàn)棣? [-50°,-36°,-23°,-12°,0°,10°,25°,40°,50°],信噪比SNR=5 dB，Monte-Carlo試驗(yàn)次數(shù)為P=30下CMR算法在不同快拍數(shù)下的RMSE的實(shí)驗(yàn)結(jié)果。由圖可知，隨著快拍數(shù)的增加，RMSE隨之減小，最后趨于平穩(wěn)，說(shuō)明算法魯棒性較好。

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

本文提出了一種基于虛擬陣列插值的矩陣重構(gòu)DOA估計(jì)算法，該算法通過(guò)陣列插值的方式將互質(zhì)陣列的虛擬陣列拓展成虛擬均勻線陣，利用均勻虛擬線陣等效信息初始化協(xié)方差矩陣，最后通過(guò)并優(yōu)化協(xié)方差矩陣的相應(yīng)矢量的首個(gè)元素，利用重構(gòu)的低秩協(xié)方差矩陣進(jìn)行DOA估計(jì)，充分利用了互質(zhì)陣列拓展的虛擬陣列的信息。將本文算法與最大連續(xù)虛擬陣元的空間平滑協(xié)方差矩陣DOA估計(jì)算法進(jìn)行對(duì)比，在不同信源數(shù)時(shí)算法的檢測(cè)空間譜、分辨率性能和不同信源數(shù)下的均方根誤差進(jìn)行了仿真實(shí)驗(yàn)，結(jié)果表明，該算法可實(shí)現(xiàn)更大的自由度和更高的分辨率，并呈現(xiàn)出較好的魯棒性。