基于ANM的寬帶非相干子空間DOA估計(jì)方法

石 娟，張群飛，毛琳琳，史文濤，王偉東

(1.西北工業(yè)大學(xué)航海學(xué)院,陜西 西安 710072；2.中國(guó)科學(xué)院聲學(xué)研究所,北京 100190)

0 引言

波達(dá)方位(Direction of Arrival，DOA)估計(jì)是陣列信號(hào)處理重要的研究?jī)?nèi)容之一，其應(yīng)用涉及雷達(dá)、通信、聲納、勘探等眾多工程領(lǐng)域[1-2]。DOA估計(jì)方法主要包括子空間分類方法和子空間擬合方法。目前，比較成熟、穩(wěn)定的DOA估計(jì)方法大都是基于陣列的空間譜估計(jì)方法，如多重信號(hào)分類方法(Multiple Signal Classification，MUSIC)和旋轉(zhuǎn)不變子空間方法(Estimation of Signal Parameters via Rotational Invariance Techniques，ESPRIT)等，以它們?yōu)榇淼目臻g譜方法突破了瑞利限，使測(cè)向定位技術(shù)實(shí)現(xiàn)了飛躍[3]。但這些方法的實(shí)現(xiàn)需要大量獨(dú)立同分布的測(cè)量數(shù)據(jù)(大快拍數(shù)據(jù))。在實(shí)際環(huán)境中，獲取大快拍數(shù)據(jù)受到時(shí)間的制約，如高速運(yùn)轉(zhuǎn)的目標(biāo)對(duì)系統(tǒng)實(shí)時(shí)性要求很高，尤其在處理超寬帶窄脈沖等短時(shí)突發(fā)數(shù)據(jù)時(shí)，接收信號(hào)經(jīng)過相干累積后只有小快拍數(shù)據(jù)乃至單快拍數(shù)據(jù)，無(wú)法直接用常規(guī)的DOA估計(jì)方法進(jìn)行方位估計(jì)。因此，研究小快拍乃至單快拍的DOA估計(jì)便是一種有效解決方案。文獻(xiàn)[4]提出了單快拍下基于匹配跟蹤的ESPRIT的方位估計(jì)方法。文獻(xiàn)[5]中提出了單快拍及多快拍下聯(lián)合Matrix Pencil 和 ESPRIT的方位估計(jì)方法。隨著壓縮感知技術(shù)的發(fā)展，稀疏重構(gòu)的思想被用于DOA估計(jì)中[6]，然而，稀疏重構(gòu)方法對(duì)網(wǎng)格劃分比較敏感，容易造成信號(hào)基失配。為了解決信號(hào)失配問題， Tang. 等人提出了基于最小原子范數(shù)(Atomic Norm Minimization，ANM)的無(wú)網(wǎng)格壓縮感知理論[7]。ANM是一種結(jié)構(gòu)優(yōu)化方法，其基于陣列流型矩陣的范德蒙結(jié)構(gòu)，通過半正定規(guī)劃(Semidefinite Programming，SDP)可以構(gòu)造恢復(fù)出所需的Toeplitz矩陣。目前，國(guó)內(nèi)外已有不少基于ANM 的DOA估計(jì)的研究，如文獻(xiàn)[8—9]，這些方法都是基于窄帶信號(hào)DOA估計(jì)的研究。

眾所周知，在實(shí)際應(yīng)用中，寬帶信號(hào)源大量存在，并且它能夠提供比窄帶更多的信息，有利于目標(biāo)信號(hào)檢測(cè)、參數(shù)估計(jì)和特征提取[10]。目前，寬帶信號(hào)處理方法的研究主要基于兩大類[10]：一類是基于非相干合成的寬帶處理方法(Incoherent Signal Subspace Method，ISM)，如文獻(xiàn)[11]提出了基于歸一化的克拉美羅界(Cramer-Rao bound, CRB)加權(quán)非相干寬帶DOA估計(jì)方法；另一類是基于相干合成的寬帶處理方法(Coherent Signal Subspace Method，CSM)，如文獻(xiàn)[12]提出了一種基于加權(quán)擬合的相干處理的子空間寬帶信號(hào)DOA估計(jì)方法。然而，這兩種寬帶處理的子空間DOA估計(jì)方法的實(shí)現(xiàn)都依賴于大快拍數(shù)據(jù)，無(wú)法直接有效的處理小塊拍數(shù)據(jù)。

本文為了解決寬帶信號(hào)在單快拍下的DOA估計(jì)失效問題，提出了基于ANM的寬帶非相干子空間DOA估計(jì)方法(Atomic Norm Minimization Incoherent Subspace Method，A-ISM)。該方法利用子帶陣列流型矩陣的范德蒙結(jié)構(gòu)，對(duì)每個(gè)子帶分別進(jìn)行ANM結(jié)構(gòu)優(yōu)化，再利用Root-Music算法來(lái)估計(jì)角度，最后將子帶估計(jì)結(jié)果的均值作為方位估計(jì)結(jié)果。

文中(·)T、(·)H分別表示轉(zhuǎn)置和共軛轉(zhuǎn)置運(yùn)算、tr(·)表示矩陣求跡運(yùn)算、‖·‖為范數(shù)運(yùn)算、‖a‖2表示a的2范數(shù)、‖A‖F(xiàn)代表矩陣A的Frobenius范數(shù)，E[·]表示期望運(yùn)算。

1 寬帶信號(hào)模型和非相干子空間方法

1.1 寬帶信號(hào)模型

假設(shè)觀測(cè)空間存在K個(gè)遠(yuǎn)場(chǎng)寬帶信號(hào)(θ1,θ2,…,θK)入射到M元均勻線陣(M>K)，其中，θk表示第個(gè)k信號(hào)的入射角，如圖1所示。假設(shè)陣列接收的信號(hào)不存在相關(guān)性，各陣元間噪聲相互獨(dú)立，且為平穩(wěn)、零均值的高斯空間白噪聲，其方差為σ2。

圖1 陣列結(jié)構(gòu)Fig.1 Array structure

文中假設(shè)信號(hào)源的個(gè)數(shù)已知,且寬帶信號(hào)的工作頻率在[fl，fh]之間，則第m個(gè)陣元在第t時(shí)刻接收到的數(shù)據(jù)可以表示為：

(1)

式(1)中，sk和wm分別表示第m陣元接收到來(lái)自第k個(gè)信號(hào)源的信號(hào)和噪聲，τm(θk)表示位于θk的目標(biāo)入射到第m陣元和相對(duì)于參考陣元的時(shí)延,L代表快拍數(shù)。

寬帶信號(hào)帶寬為B=fh-fl，將其劃分為N個(gè)子頻帶(窄帶),各子頻帶的頻率分別為f0,f1,…,fN-1，則第n個(gè)子帶的接收數(shù)據(jù)為：

Yn=A(fn,θ)Sn+Wn
n=1,2,…,N

(2)

式(2)中，Sn=[s1s2…sK]T∈K×L，Wn∈M×L。為了簡(jiǎn)化表示，在下面的章節(jié)中A(fn,θ)用An(θ)表示，且An(θ)∈M×K是范德蒙矩陣。第n個(gè)子帶的導(dǎo)向矢量矩陣表示為：

An(θ)=[a(fn,θ1)a(fn,θ2) …a(fn,θK)]

(3)

其中，第k個(gè)信號(hào)在第n頻帶上的導(dǎo)向矢量為：

(4)

式(4)中，d表示陣列陣元間距，c為信號(hào)的傳播速度。因此，寬帶信號(hào)的的模型可以寫成：

(5)

1.2 非相干子空間方法

基于非相干信號(hào)ISM方法是最早出現(xiàn)的寬帶DOA估計(jì)方法，主要思路是將寬帶信號(hào)劃分為多個(gè)子帶，對(duì)每個(gè)子帶利用信號(hào)子空間或噪聲子空間， 求出空間譜或真實(shí)信號(hào)的根來(lái)估計(jì)方位角度，最后將各子帶所得結(jié)果平均作為最終方位估計(jì)結(jié)果。

根據(jù)寬帶信號(hào)的模型，第n個(gè)子帶的協(xié)方差矩陣為：

(6)

(7)

式(7)中，ynl表示第n子帶的第l次快拍數(shù)據(jù)。當(dāng)L

圖2 平滑方法原理圖Fig.2 Smoothing method principle structure

M元均勻線陣分成相互交錯(cuò)的P個(gè)子陣(P=1+M/2)，每個(gè)子陣的個(gè)數(shù)為M/2，第n子帶平滑后的數(shù)據(jù)可以表示為：

(8)

平滑后協(xié)方差矩陣表示為：

(9)

(10)

其中，特征值λi?λg，vi和vg分別代表特征值λi和λg對(duì)應(yīng)特征向量，且信號(hào)子空間Us和噪聲子空間Uw可描述為：

(11)

在得到噪聲子空間Uw后，將其應(yīng)用到Root-Music方法中進(jìn)行DOA估計(jì)。首先，Root-Music方法定義為如下一個(gè)多項(xiàng)式[13]：

(12)

其中

Q(z)=[1z…zM-1]T

(13)

由式(12)求出K個(gè)接近單位圓的根z1,z2,…,zK。

然后, 根據(jù)下式得到第n個(gè)子帶的DOA估計(jì)角度

(14)

隨后，對(duì)N個(gè)子帶估計(jì)出來(lái)的角度求均值，得到最終DOA估計(jì)：

(15)

2 基于ANM的寬帶非相干子空間方位估計(jì)方法

ISM方法在處理寬帶信號(hào)時(shí)，應(yīng)用了窄帶平均思想，并且要求各子帶數(shù)據(jù)的快拍數(shù)必須大于陣元個(gè)數(shù)，即L>M，否則不能滿足子空間類方法對(duì)空間協(xié)方差矩陣滿秩的要求。當(dāng)數(shù)據(jù)有限，尤其快拍數(shù)L=1時(shí)，直接應(yīng)用傳統(tǒng)的ISM方法進(jìn)行DOA估計(jì)時(shí)，系統(tǒng)完全不能工作。雖然平滑ISM方法可以估計(jì)方位，但是必須以犧牲目標(biāo)分辨率為代價(jià)[14]。針對(duì)單快拍下的方位估計(jì)，本文提出了基于最小原子范數(shù)的寬帶非相干子空間DOA估計(jì)方法(A-ISM)。該方法可以直接用單快拍數(shù)據(jù)，通過ANM構(gòu)造并會(huì)恢復(fù)出Toeplitz 矩陣進(jìn)行DOA估計(jì)。

在單快拍L=1條件下，第n子帶接收的數(shù)據(jù)可表示為：

yn=xn+wn=
An(θ)sn+wnn=1,2,…，N

(16)

式(16)中，yn為M×1的向量，xn∈M×1和wn∈M×1分別表示信號(hào)向量和噪聲向量。首先考慮無(wú)噪聲條件，即式(16)中只包含了與角度有關(guān)的信號(hào)信息，其可寫成：

yn=xn=An(θ)snn=1,2,…,N

(17)

原子核定義為：

(18)

式(17)表明，xn是An(θ)的線性組合，因此，xn的原子范數(shù)可表示成：

(19)

根據(jù)文獻(xiàn)和式(4)可知，如果信號(hào)源是稀疏的，則xn中的sin(θ)就是稀疏的，式(17)就能夠進(jìn)行稀疏原子分解。求解原子范數(shù)最小化是一個(gè)NP(Non-Deterministic Polynormial)問題，很難直接求解。而原子范數(shù)具有半正定規(guī)劃性質(zhì)(SDP)，采取線性半正定規(guī)劃方法就能在多項(xiàng)式時(shí)間內(nèi)求解原子范數(shù)；同時(shí)，由Caratheodory-Toeplitz引理可知，任何半正定Toeplitz矩陣都能進(jìn)行范德蒙分解[8]。因此，可以將原子范數(shù)最小化問題轉(zhuǎn)化為半正定規(guī)劃問題，其可以描述為[8]：

(20)

由式(20)可以解出一個(gè)最優(yōu)的向量un=[un1un2…unM]T，用以構(gòu)造最優(yōu)的Toeplitz矩陣[7]

(21)

式(20)中的半正定的約束，根據(jù)舒爾補(bǔ)(Schur Complement)條件，其可等價(jià)于[14]：

(22)

(23)

實(shí)際應(yīng)用中，噪聲總是伴隨著信號(hào)存在，因此，式(20)中必須考慮對(duì)噪聲作用進(jìn)行約束。噪聲存在時(shí)，ANM的SDP可以重新描述為：

(24)

式(24)中，ρ是一個(gè)正則化參數(shù)，用來(lái)調(diào)節(jié)T(un)的Toeplitz結(jié)構(gòu)與yn中噪聲方差。給定了噪聲方差σ2，ρ可以根據(jù)下式計(jì)算[6]：

(25)

3 仿真與分析

本章以水下被動(dòng)探測(cè)系統(tǒng)為背景，在單快拍下，利用計(jì)算機(jī)仿真分別從分辨概率、方位估計(jì)歸一化均方誤差對(duì)A-ISM方法和平滑ISM方法的DOA估計(jì)性能進(jìn)行比較分析。

仿真實(shí)驗(yàn)中均采用16元均勻線陣，目標(biāo)源個(gè)數(shù)K=2，工作頻率為3 000～3 100 Hz，子頻帶個(gè)數(shù)N=3，快拍數(shù)L=1。陣元分布均滿足陣元間距為寬帶信號(hào)最低頻率fl對(duì)應(yīng)波長(zhǎng)的半波長(zhǎng)，水聲傳播速度c=1 500 m/s。仿真中利用CVX包解決式(20)和式(24)中的SDP問題。

首先，圖3比較了目標(biāo)源位于(30° 36°)時(shí)，本文所提的A-ISM方法和平滑ISM方法的分辨概率。本次仿真實(shí)驗(yàn)的分辨概率是指在若干次蒙特卡羅實(shí)驗(yàn)中，能夠正確分辨兩個(gè)或多個(gè)目標(biāo)的概率。正確分辨需滿足以下條件：

(26)

式(26)中，B=arcsin(0.44λ/Md) 。由圖3可以得出，分辨概率為95%時(shí)，A-ISM方法需要的信噪比為4 dB，而平滑ISM方法需要的信噪比為12 dB，A-ISM方法相較平滑ISM方法，信噪比提高了8 dB。

圖3 分辨概率與信噪比的關(guān)系，L=1,θ=(30° 36°)Fig.3 Probability of resolution versus SNRs，L=1， θ=(30° 36°)

圖4中依然假設(shè)目標(biāo)源位于(30° 36°)，給出了A-ISM方法和平滑ISM的方位估計(jì)性能隨信噪比的變化情況，且本節(jié)的仿真實(shí)驗(yàn)中，方位估計(jì)選取歸一化均方誤差(nMSE)作為衡量方位估計(jì)誤差的標(biāo)準(zhǔn)，定義為：

(27)

圖4 歸一化均方根誤差與信噪比的關(guān)系，L=1，θ=(30° 36°)Fig.4 nMSE of DOA versus SNRs，L=1， θ=(30° 36°)

圖5和圖6比較了兩種方法在不同信噪比下，歸一化均方誤差隨著兩目標(biāo)夾角Δθ的變化情況，以此來(lái)衡量方法的方位估計(jì)性能。從圖中可以看出，平滑ISM方法與A-ISM方法的歸一化均方誤差曲線存在一定的差距，后者估計(jì)性能更優(yōu)。再比較圖5和圖6，隨著信噪比的提高，兩種方法的歸一化均方誤差均有減小，但相較平滑ISM方法，A-ISM方法的歸一化均方誤差依舊較小，表現(xiàn)出較優(yōu)的估計(jì)性能。

圖5 歸一化均方誤差與的目標(biāo)夾角關(guān)系,SNR=10 dB,L=1Fig.5 nMSE of DOA versus the DOA separation,L=1

圖6 歸一化均方誤差與的目標(biāo)夾角關(guān)系, SNR=15 dB,L=1Fig.6 nMSE of DOA versus the DOA separation，SNR=15dB，L=1

4 結(jié)論

本文提出了基于ANM的寬帶非相干子空間DOA估計(jì)方法(A-ISM)。該方法通過各子帶進(jìn)行ANM優(yōu)化處理，恢復(fù)出最優(yōu)的Toeplitz矩陣，該矩陣經(jīng)過特征分解能準(zhǔn)確區(qū)分信號(hào)子空間和噪聲子空間，再結(jié)合Root-Music算法，最終得到誤差較小的估計(jì)角度。該方法不僅解決了單快拍下經(jīng)典寬帶信號(hào)DOA估計(jì)方法失效問題，而且提高了方位估計(jì)性能。仿真結(jié)果驗(yàn)證了所提DOA估計(jì)方法的正確性和有效性。