羅雪

(北京理工大學(xué) 信息與電子學(xué)院，北京 100081)

多信源測向在雷達(dá)、聲納、移動通信等領(lǐng)域受到了極大關(guān)注. 以多重信號分類(multiple signal classification, MUSIC)[1]和基于旋轉(zhuǎn)不變技術(shù)的信號參數(shù)估計(jì)(estimating signal parameter via rotational invariance techniques,ESPRIT)[2]為代表的超分辨算法，能夠利用子空間理論識別入射角度相距很近的信源，因此在過去幾十年間被廣泛應(yīng)用[3]. 盡管特征結(jié)構(gòu)方法有很多潛在優(yōu)勢，但由于它們對陣列流形信息的精度有極高要求，將其直接應(yīng)用在實(shí)際系統(tǒng)中是很困難的. 換句話說，超分辨算法的性能嚴(yán)重依賴于陣列流形的精確性. 實(shí)際上，陣列流形不可避免地會受到互耦、位置變動和陣列增益或相位不確定性的影響，這將造成陣列接收信號明顯的幅度與相位失真. 直接用基于特征結(jié)構(gòu)的方法將會造成測向性能嚴(yán)重惡化[4-7]. 陣列校準(zhǔn)是一個能有效緩解波達(dá)方向(direction of arrival, DOA)估計(jì)偏差的方法.

事實(shí)上，在過去幾十年間，人們已經(jīng)研究了多種校準(zhǔn)方法. 包括插值法[8-9]，已知位置校準(zhǔn)信源的最大似然法[10]，最小均方法[11]，自校準(zhǔn)方法[12-14]以及最大后驗(yàn)估計(jì)器[15-16]. 在近些年，多個基于互耦矩陣(mutual coupling matrix, MCM)特征的算法被發(fā)展. YE[17]通過在陣列兩側(cè)設(shè)置傳感器獲得了DOA的初始估計(jì). 因此，運(yùn)用原始陣列數(shù)據(jù)能夠直接執(zhí)行空間譜的一維搜索. 此外，還可以將估計(jì)的互耦參數(shù)用于補(bǔ)償，搜索結(jié)果能夠通過迭代進(jìn)一步細(xì)化. 從中心子陣收到的信號被直接運(yùn)用到傳統(tǒng)的MUSIC，因而MCM被假設(shè)有一個復(fù)對稱的Toeplitz矩陣[18]. 為了提升DOA估計(jì)的性能，LIAO等[19]中提出利用MCM的特殊結(jié)構(gòu)參數(shù)化了方向向量，從而用整個陣列來實(shí)現(xiàn)DOA估計(jì)，并通過互耦補(bǔ)償提升了結(jié)果的精確性. 此外，這些校準(zhǔn)方法也被應(yīng)用于遠(yuǎn)場信源[20]和非圓信號的DOA問題中[21].

本文假設(shè)MCM是有有限個非零元素的帶對稱Toeplitz矩陣. 通過運(yùn)用文獻(xiàn)[19]中參數(shù)化的方向向量，能夠用整個陣列獲得初步的DOA估計(jì). 之后，估計(jì)角度依賴系數(shù)(angularly dependent coefficients , ADCs)，利用子空間的正交性和有全部ADCs信息的互耦系數(shù)進(jìn)行進(jìn)一步的互耦補(bǔ)償.

1 信號模型與問題分析

本文考慮了由M個傳感器組成的均勻線陣(uniform linear array, ULA)收到了K個窄帶遠(yuǎn)場信源s1(t),s2(t),…,sK(t). 信源分別從與陣列法線相隔θ1,θ2,…,θK角度的方向入射. 假設(shè)相鄰傳感器間的距離為d，載波波長為λ，由于傳感器間的相互作用，互耦影響不能夠被忽略.M×1維的陣列輸出向量可以表示為

x(t)=CAs(t)+n(t)

(1)

式中：x(t)=[x1(t)x2(t)…xM(t)]T；s(t)=[s1(t)s2(t)…sK(t)]T；n(t)=[n1(t)n2(t)…nM(t)]T分別表示了接收信號向量，源信號向量和噪聲向量;符號[·]T表示轉(zhuǎn)置;A=[a(θ1)a(θ2)…a(θK)]是陣列流形矩陣，其中的變量a(θk)=[1β(θk)β2(θk)…βM-1(θk)]T是方向向量，β(θk)=exp{-j2πdsinθk/λ};C∈M×M是MCM，通常被認(rèn)為獨(dú)立于信號DOA[22]，代表了陣列中任意兩個傳感器間的相互作用. 假設(shè)加性傳感器噪聲ni(t),i=1,2,…,M為獨(dú)立同分布的高斯白噪聲，有相同的方差σ2.si(t)和ni(t)為零均值廣義平穩(wěn)隨機(jī)過程.

關(guān)于耦合模型的幾個研究[12,23]表明一對傳感器間的耦合幾乎是相同的. 因此，ULA的MCM是一個帶狀對稱的Toeplitz矩陣. 此外，兩個傳感器間的互耦與它們的距離成反比，如果兩個傳感器相距幾個波長，系數(shù)將變?yōu)?. 用cij=cji=c|i-j|表示ULA的第i個陣元與第j個陣元間的互耦系數(shù)，并假設(shè)MCM中共有P個不同的非零元素，其中自耦c0被歸一化為1. 那么，M×M維的矩陣C能夠被表示為

C=Toeplitz{1,c1,…,cP-1,0M-P}

(2)

其中，Toeplitz{·}表示了一個對稱的Toeplitz矩陣，由P×1維向量c=[1c1…cP-1]T構(gòu)造. 定義方向θ的等效方向向量為

am(θ)=Ca(θ)

(3)

從公式(1)中可見，接收信號的協(xié)方差矩陣為

Rx=E[x(t)xH(t)]=CARsAHCH+σ2I=

(4)

式中：(·)H表示了共軛轉(zhuǎn)置操作；Rs=E[s(t)sH(t)]為s(t)的協(xié)方差矩陣；Am=[am(θ1)am(θ2)…am(θk)]是等效的陣列流形矩陣. 如文獻(xiàn)[2]中所證明的，運(yùn)用Rx的特征值分解(eigenvalue decomposition, EVD)，式(4)能夠被重新寫為

Rx=[Es,En]Λ[Es,En]H

(5)

其中Λ=diag{Λs}為一個矩陣，M個特征值降序排列在對角線上，其余位置為0.Es∈M×K和En∈M×(M-K)分別由K個最大特征值與M-K個較小特征值對應(yīng)的特征向量組成. 根據(jù)著名的子空間的算法MUSIC[1]，Am=CA的列張成的信號子空間與En張成的噪聲子空間是正交的. 因此下式成立

(6)

2 互耦補(bǔ)償與波達(dá)方向估計(jì)

首先通過參數(shù)化MCM重新制定了等效的方向向量. 之后，在存在未知互耦的情況下，通過應(yīng)用整個陣列采樣數(shù)據(jù)解決了DOA估計(jì)問題，獲得更好的性能[19,24]. 最后，為了進(jìn)一步細(xì)化，本文提出了一種新的互耦補(bǔ)償算法，該算法利用了MCM的特殊結(jié)構(gòu)，沒有信息損失.

2.1 利用全部陣列采樣的DOA估計(jì)

定義rk為M×1維等效方向向量am(θ)的第k個元素. 結(jié)合公式(2)(3)，am(θ)能夠被表示為

am(θ)=[r1r2…rM]T

(7)

其中

am(θ)=g0Γ(θ)a(θ)

(8)

其中，Γ(θ)=diag{μ1,…,μP-1,1,…,α1,…,αP-1}，中間有M-2(P-1)個1，且

μk=rkβ1-k(θ)/g0,k=1,2,…,P-1

(9)

αk=rk+(M-P+1)β1-k(θ)/g0,k=1,2,…,P-1

(10)

是ADCs，它們?yōu)閏k和β(θ)的函數(shù). 式(8)～(10)表明，與角度無關(guān)的MCM被轉(zhuǎn)換成角度有關(guān)的表達(dá)式g0Γ(θ). 因?yàn)棣?θ)是有2(P-1)個未知變量的對角矩陣，a(θ)為一個向量，交換他們的元素并將Γ(θ)中1的數(shù)量減少為1個[19]，可得

am(θ)=g0T(θ)v(θ)

(11)

v(θ)=[μ1…μP-11α1…αP-1]T為一個(2P-1)×1維向量，并且

(12)

是一個M×(2P-1)維矩陣，由一個(M-2P+2)×1維向量和兩個(P-1)×(P-1)維對角矩陣組成.

考慮式(6)中的子空間特性和式(11)中am(θ)的定義，能夠在g0≠0的條件下得到如下的等式：

vH(θ)Q(θ)v(θ)=0

(13)

(14)

要注意的是，Q(θ)并不包含c的任何信息. 因此即使在互耦未知的情況下也能使用這個譜函數(shù). 與Svantesson提出的算法[17]相比，這個譜函數(shù)利用了全部的陣列采樣，不需要提取中間陣列進(jìn)行互耦消除. 因此譜估計(jì)沒有任何信息損失. 此方法僅適用于如下條件被滿足的情況：

K≤M-2P+1

(15)

在這個情況下，(M-K)×(2P-1)維矩陣EHT(θ)是滿列秩的.

2.2 互耦補(bǔ)償

(16)

(17)

Q′(∶,2:2P-1)v′(2∶P)=-qP

(18)

因此，通過求解上述方程能夠得到μk和αk的估計(jì)：

v′(2∶P)=-Q′(∶,2∶2P-1)?qP

(19)

(·)?代表了偽逆操作.

vg=

(20)

那么式(21)成立：

(21)

其中：

(22)

(23)

結(jié)合式(20)、(23)，互耦參數(shù)向量c′能夠通過求解式(21)得到

(24)

(25)

通過最小二乘法求解公式(25)，可以得到c′的精確估計(jì)為

(26)

上面的方法提供了一個互耦補(bǔ)償?shù)姆绞?，即?dāng)向量c′被確定后，矩陣C可以通過將其元素定位在相應(yīng)的對角線上形成. 因此，DOA估計(jì)能夠被獲得通過搜索譜函數(shù)的峰值：

(27)

DOA估計(jì)性能可以被進(jìn)一步的提升通過重復(fù)上述過程. 文中提出的DOA估計(jì)與互耦補(bǔ)償算法總結(jié)如下.

(1) 在t=t1,t2,…,tL時刻陣列收到信號x(t)的L個快拍，得到數(shù)據(jù)矩陣：

X=[x(t1)x(t2)…x(tL)]

(28)

(2) 運(yùn)用上述的數(shù)據(jù)矩陣生成協(xié)方差矩陣：

Rx=XXH/L

(29)

(3) 對Rx進(jìn)行特征值分解，得到噪聲子空間En.

(7) 用估計(jì)的C和式(27)提升DOA估計(jì)性能.

(8) 重復(fù)第(5)～(7)步可以得到更精確的方向估計(jì).

3 波達(dá)方向估計(jì)性能對比

考慮兩個獨(dú)立信源從遠(yuǎn)場入射到ULA，有入射角度θ1=-10°和θ2=20°. 陣列中的傳感器間隔為d=λ/2. 有效互耦參數(shù)的數(shù)量為P=3. 設(shè)置c=[10.433 01-0.25i0.141 42-0.141 42i]，它被用在了YE[17]中的第2個仿真，保證了在任意方向θ有g(shù)0≠0.

3.1 無互耦補(bǔ)償下的估計(jì)性能比較

第一個實(shí)驗(yàn)評估了沒有互耦補(bǔ)償?shù)那闆r下第(4)步中的DOA估計(jì)性能. 假設(shè)陣列中陣元數(shù)為M=7，快拍數(shù)為500. 均方根誤差(RMSE)被用于比較不同算法的DOA估計(jì)精度，它通常能夠被計(jì)算為

(30)

RMSE作為信噪比的函數(shù)如圖1(a)所示. 在低信噪比時，由于用到了整個陣列，本文提出的算法優(yōu)于文獻(xiàn)[7]中的方法. 隨著信噪比增加，所有方法DOA估計(jì)的RMSE都逐步降低. 當(dāng)信噪比高于10 dB時，用于矩陣C未知情況下的兩種方法有幾乎相同的估計(jì)準(zhǔn)確性. 圖1(b)顯示了陣列尺寸對于RMSE的影響，信源數(shù)M的選擇應(yīng)該滿足式(15). 能夠看到在M≤10時本文的方法稍優(yōu)于文獻(xiàn) [17]中的方法.

圖1 DOA估計(jì)的RMSE作為SNR的函數(shù)和DOA估計(jì)的RMSE作為陣元數(shù)M的函數(shù)Fig.1 RMSE of DOA versus SNR and RMSE of DOA versus the array number M

運(yùn)用初始DOA估計(jì)可以得到互耦向量c. 在第二個仿真中，考慮與第一個仿真相同的情況，定義c的RMSE為

(31)

圖2 互耦矩陣(MCM)估計(jì)的RMSE作為SNR和快拍數(shù)的函數(shù)Fig.2 RMSE of mutual coupling matrix(MCM) estimates and RMSE of MCM estimates versus snapshots

3.2 有互耦補(bǔ)償下的估計(jì)性能比較

考慮兩個信源分別從角度θ1和θ2=θ1+Δθ入射到陣元數(shù)為M=7的ULA上，RSN=0 dB，快拍數(shù)L=500. 定義：

|θ1-θ1|+|θ2-θ2|<|θ2-θ1|

圖3 DOA的可識別率作為角度間隔的函數(shù)以及DOA估計(jì)的RMSE分別作為SNR、快拍數(shù)、陣列尺寸的函數(shù)Fig.3 Probability of DOA identification versus the angle interval RMSE of DOA versus SNR,snapshots and the array number

4 結(jié) 論

本文解決了存在未知互耦時的DOA估計(jì)問題. 該方法基于子空間理論，在沒有校準(zhǔn)源和輔助傳感器的情況下首先用整個陣列估計(jì)初始DOA. 假設(shè)在空間中不存在盲角，該方法通過從角度相關(guān)系數(shù)間接估計(jì)耦合系數(shù)來進(jìn)行互耦補(bǔ)償. 最后，在不丟棄任何ADCs的情況下，通過求解最小二乘問題來確定互耦系數(shù). 仿真結(jié)果表明，提出的算法在低SNR和小尺寸陣列的情況下可獲得更好的性能. 同時，提出方法的魯棒性也得到了驗(yàn)證.