虞 飛，宋 俊，余 赟，蘇 冰

(海軍研究院，北京 100071)

0 引 言

基于稀疏重構(gòu)理論的陣列測向方法是近十年陣列信號處理領(lǐng)域新發(fā)展起來的十分活躍的一個研究分支，已取得了豐碩的理論成果[1-4]。這類方法在傳統(tǒng)的基于子空間的超分辨算法的基礎(chǔ)上，引入了目標(biāo)信號的空間稀疏信息，從而大大提高了目標(biāo)信號的到達(dá)角估計(jì)性能。傳統(tǒng)的子空間類到達(dá)角估計(jì)算法要求在高信噪比、非相關(guān)、大樣本、理想陣列流形條件下才可獲得優(yōu)異的性能[5-7]，盡管近年來圍繞其中某一種非理想因素進(jìn)行算法優(yōu)化的研究成果很多[8-10]，但是算法的實(shí)用性因?qū)嶋H使用場景下多個非理想因素往往同時存在而大大受限。

稀疏重構(gòu)類到達(dá)角估計(jì)方法之所以能夠成為近年來的研究熱點(diǎn)，關(guān)鍵原因在于此類方法在低信噪比、相干背景、小樣本、陣列未校準(zhǔn)等不利條件下都具有良好的目標(biāo)信號測向精度和多目標(biāo)分辨能力，相對于傳統(tǒng)的子空間類算法在實(shí)用性方面具有更大潛力[11]。

雖然稀疏重構(gòu)類到達(dá)角估計(jì)方法在整體魯棒性方面優(yōu)于子空間類估計(jì)技術(shù)，但是前者的主流算法在計(jì)算復(fù)雜度方面處于相對劣勢，而且還涉及額外的正則化參數(shù)選取問題。本文通過稀疏重構(gòu)得到傳感器陣列輸出數(shù)據(jù)的稀疏表示模型，提出了一種基于最小均方誤差(Minimum Mean-Square Error,MMSE)準(zhǔn)則迭代實(shí)現(xiàn)的單快拍到達(dá)角估計(jì)算法(Iterative Implementation of MMSE, II-MMSE)。該算法采用MMSE框架，考慮了陣列接收數(shù)據(jù)中的模型“噪聲”協(xié)方差信息，對于陣列模型誤差具有較好的魯棒性。II-MMSE算法保持了稀疏重構(gòu)類到達(dá)角估計(jì)方法在諸多非理想因素下的整體魯棒性優(yōu)勢，克服了主流稀疏重構(gòu)類方法的高計(jì)算復(fù)雜度問題，而且避免了正則化參數(shù)選取，進(jìn)一步推動了稀疏重構(gòu)類到達(dá)角估計(jì)方法的工程實(shí)用化。

1 陣列輸出模型的稀疏重構(gòu)

考慮K個遠(yuǎn)場窄帶平面波入射到由M個無指向性單元構(gòu)成的傳感器陣列，且K個信號與傳感器陣列位于同一平面內(nèi)。將目標(biāo)信號可能到達(dá)角范圍空間Θ進(jìn)行N次均勻網(wǎng)格劃分，可將傳統(tǒng)的陣列輸出模型轉(zhuǎn)化為稀疏重構(gòu)模型[11]：

2 基于MMSE準(zhǔn)則迭代實(shí)現(xiàn)的單快拍DOA估計(jì)算法

根據(jù)最小均方誤差(MMSE)準(zhǔn)則和式(1)中的稀疏重構(gòu)模型，取目標(biāo)函數(shù)f(t)：

式中：W(t)為M×N維復(fù)權(quán)矩陣，γ(t)顯然是經(jīng)過M×N型濾波器組濾波后的期望輸出信號矢量。根據(jù)Wiener濾波器理論相關(guān)結(jié)論，容易求得式(2)中的目標(biāo)函數(shù)在MMSE意義上的最優(yōu)權(quán)矩陣為

假定γ(t)與n(t)之間相互統(tǒng)計(jì)獨(dú)立，則由式(1)可得：

式中：Rn=E[n(t)nH(t)]為M×M維噪聲協(xié)方差矩陣，對于功率為σ2n的空間高斯白噪聲，Rn可以簡化為Rn=σ2nIM，P=E[γ(t)γH(t)]為N×N維稀疏信號協(xié)方差矩陣。假設(shè)稀疏信號矢量γ(t)內(nèi)各個分量互不相關(guān)，則有：

其中：pn=E[γn(t)γ*n(t)]，n=1，…，N，故P可稱為稀疏功率對角陣。

由式(4)和(5)，可將式(3)進(jìn)一步表示為

采用空間匹配濾波器(Matched Filter, MF)即常規(guī)波束形成(Conventional Beamforming, CBF)算法，可以得到稀疏信號γ(t)的粗略估計(jì)：

(t)可作為MMSE迭代算法中稀疏信號γ(t)的初始估計(jì)值：

則稀疏功率對角陣的初始值為

式中：⊙表示矩陣的哈達(dá)瑪(Hadamard)乘積，即矩陣對應(yīng)元素乘積。再由式(7)可得MMSE迭代算法中權(quán)矩陣的更新值：

式中：下標(biāo)i表示第i次迭代。稀疏信號γ(t)的最小均方誤差估計(jì)值為

類似于式(10)，可得第i次迭代時稀疏功率對角陣的估計(jì)為

式(11)、(12)和(13)構(gòu)成了MMSE迭代算法第i次迭代的實(shí)現(xiàn)過程。當(dāng)某一次迭代估計(jì)值滿足ε為某一預(yù)先指定的較小正數(shù)時，算法停止迭代，或者當(dāng)算法達(dá)到預(yù)設(shè)的最大迭代次數(shù)時也停止迭代。此時，由稀疏功率對角陣的估值構(gòu)成的列向量diag[(t)]關(guān)于到達(dá)角網(wǎng)格的稀疏功率譜中，第K個譜峰對應(yīng)的角度即為目標(biāo)信號到達(dá)角估計(jì)值

綜合以上分析，現(xiàn)將MMSE迭代算法的完整流程歸納如下：

(1) 初始化

對于t時刻的陣列接收數(shù)據(jù)x(t)，有

(2) 迭代過程

(3) 算法結(jié)果

在列向量diag[(t)]關(guān)于到達(dá)角網(wǎng)格的稀疏功率譜中，第k個譜峰對應(yīng)的角度值即為目標(biāo)信號到達(dá)角估計(jì)值

3 算法收斂性分析

下面分析MMSE迭代算法的收斂性。為了說明MMSE算法迭代的收斂性，首先建立稀疏信號矢量估計(jì)的遞推表達(dá)式。由式(11)、(12)和(13)可得：

式中：“(·)?”表示矩陣的穆爾-彭羅斯(Moore-Penrose)逆。在實(shí)際應(yīng)用中，噪聲協(xié)方差矩陣Rn不能忽略。迭代過程式(15)實(shí)際上是如下加權(quán)最小范數(shù)優(yōu)化問題的一種遞推形式[12]：

式中：W為N×N維權(quán)矩陣，且W為對角陣，其迭代形式為Wi=diag[(t)]。在每一步迭代中，目標(biāo)函數(shù)都有如下關(guān)系：

其中：wk為權(quán)矩陣W的第k個對角元素，γk(t)為稀疏信號矢量γ(t)的第k個元素。

由式(17)可以看出，權(quán)矩陣W中某個對角元素相對越大，γ(t)中的相應(yīng)位置元素對目標(biāo)函數(shù)最小化的貢獻(xiàn)就相對越小，即懲罰值越小。反之亦然。因此，如果稀疏基矩陣Φ中的某一列相對于其他列來說能更好地匹配陣列測量數(shù)據(jù)x(t)，那么γi-1(t)中的對應(yīng)位置元素迭代到下一步時將得到更大值。這樣，通過設(shè)定一個可行的初始化稀疏信號矢量估計(jì)，如可使目標(biāo)函數(shù)(17)在最小化的迭代過程中，逐漸強(qiáng)化γ(t)中某一個大小相對突出的元素，同時逐漸抑制剩余元素的大小，直至γ(t)達(dá)到預(yù)設(shè)的估計(jì)精度或者這些受抑制元素近似全為0，此時算法收斂，停止迭代。這時，僅選取稀疏基矩陣Φ中的某一列來最佳地匹配陣列測量數(shù)據(jù)x(t)。

由式(17)容易得出第i次迭代時的目標(biāo)函數(shù)為

綜上分析可知，式(18)中的遞推問題收斂到最稀疏解意味著當(dāng)i→+∞且(t)≠0時，有而MMSE迭代算法的收斂條件可以進(jìn)一步寫成：

因此，如果當(dāng)i→+∞且(t)≠0時，有成立，則也成立，從而說明MMSE迭代算法是收斂的。

應(yīng)注意，MMSE迭代算法并非在任意初始化條件下的收斂都是有意義的，例如當(dāng)(t)=0時，γ(t)每步迭代的結(jié)果都為0。因此，不失一般性，可假定初始化稀疏信號矢量(t)的非0元素個數(shù)始終為N。另外，權(quán)矩陣W的對角元素wk=0意味著，通過運(yùn)算ΦWi使稀疏基矩陣Φ中的相應(yīng)列變成了0向量，說明相應(yīng)的子空間被排除出了信號子空間。

4 陣列模型誤差下的算法性能分析

在實(shí)際應(yīng)用中，由于傳感器基陣工作環(huán)境中的介質(zhì)擾動、陣列長期未校準(zhǔn)或存在校準(zhǔn)誤差、陣列有限采樣引起的幅相量化誤差等，都可能導(dǎo)致傳感器基陣模型產(chǎn)生誤差。結(jié)合式(1)，含有陣列模型誤差的陣列輸出響應(yīng)一般可以定義為

式中：z為M×1維未知的陣列模型誤差矢量，其第m個元素可表示為

其中：Δam和Δφm分別為陣元m上的隨機(jī)幅度誤差和隨機(jī)相位誤差。假設(shè)各陣元具有獨(dú)立且同分布的零均值隨機(jī)幅值誤差和隨機(jī)相位誤差，則各陣元之間的模型誤差是互不相關(guān)的，但具有相同的模型誤差的方差σ2z。式(20)還可以表示為

式中：nz(t)=[Φγ(t)]⊙(z-IM×1)等效為陣列模型誤差引起的“噪聲”矢量，這里IM×1為M×1維全1向量。由模型誤差的假設(shè)可知nz(t)為M×1維零均值矢量。

類似于式(11)，可得存在陣列模型誤差時，MMSE迭代算法的權(quán)矩陣更新為

式中：Rnz=E[nz(t)nHz(t)]表示模型噪聲協(xié)方差矩陣，且有：

其中：Z=diag[z1，…，zM]-IM。在式(24)的推導(dǎo)過程中，還利用到了Hadamard乘積的運(yùn)算性質(zhì)A⊙B=B⊙A和(A⊙B)H=AH⊙BH[13]。根據(jù)式(24)的結(jié)果，則權(quán)矩陣的更新變?yōu)?/p>

在式(25)中，Rn為僅依賴于噪聲的固定協(xié)方差矩陣，而σ2zIM⊙[H]是依賴于信號功率更新估計(jì)的自適應(yīng)噪聲協(xié)方差矩陣。當(dāng)有高功率信號源存在時，可能導(dǎo)致算法對噪聲功率的低估，或者存在陣列模型誤差等情形時，都可能會引起小的偽峰，而模型噪聲協(xié)方差項(xiàng)σ2zIM⊙[H]的存在為信號源的功率估計(jì)建立了一個可接受的動態(tài)范圍，恰好消除了這些偽峰的影響，從而使本文算法具有自適應(yīng)能力，表明MMSE迭代算法對陣列模型誤差具有較好的魯棒性。

為了增強(qiáng)算法的自適應(yīng)能力，可對式(25)進(jìn)行以下修正[14]：

其中，通過引入尺度因子α(0＜α＜1)降低了噪聲協(xié)方差項(xiàng)Rn，可以進(jìn)一步增強(qiáng)算法的自適應(yīng)能力，使自適應(yīng)噪聲協(xié)方差項(xiàng)在權(quán)矩陣更新過程中的相對貢獻(xiàn)更大。

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

考慮兩個等功率的相干窄帶遠(yuǎn)場平面波分別從不同方向入射到由12個傳感器陣元按照半波長間距布陣構(gòu)成的均勻線列陣，陣列對空間信號進(jìn)行單快拍采樣。定義信噪比RSN=10lg(Psσ2n)，如無特殊說明，實(shí)驗(yàn)中取RSN=10 dB。

實(shí)驗(yàn)中，第l個陣元對應(yīng)的模型誤差可表示為

式中：ρ為相對于標(biāo)準(zhǔn)差的百分比，N(0，1)為零均值，單位方差的實(shí)高斯隨機(jī)分布噪聲。由式(21)可知，E[zl]=1，則故可以通過zl-1的1 000次獨(dú)立實(shí)現(xiàn)來估計(jì)zl的方差，即。顯然，ρ越大，模型誤差的方差也越大，而ρ=0、α=1時，算法退化為理想陣列模型。如無明確說明，實(shí)驗(yàn)中噪聲協(xié)方差尺度因子α=1/8。

5.1 算法收斂性實(shí)驗(yàn)

假設(shè)兩個目標(biāo)信號分別以到達(dá)角參數(shù)θ1=-10°，θ2=60°入射到上述傳感器陣列，陣列模型誤差百分比ρ=10。圖1為II-MMSE算法分別在初始化、迭代1次、5次和10次過程中對目標(biāo)信號到達(dá)角的歸一化稀疏功率譜圖。這里算法初始化結(jié)果為通過空間匹配濾波器估計(jì)所得。

圖1 II-MMSE算法迭代不同次數(shù)得到的兩個信號到達(dá)角估計(jì)的歸一化功率譜Fig.1 Normalized power spectrums of the DOA estimate for two signals obtained by the II-MMSE algorithm with different times of iterations

從仿真結(jié)果可發(fā)現(xiàn)，盡管存在陣列模型誤差，II-MMSE算法仍能在真實(shí)目標(biāo)方向角上形成尖銳譜峰，表明算法對目標(biāo)信號具有良好的到達(dá)角估計(jì)精度和多目標(biāo)分辨能力。當(dāng)算法迭代到一定次數(shù)時(一般不超過15次)，功率譜圖的偽峰及旁瓣逐漸被抑制，且目標(biāo)方向?qū)?yīng)的譜峰十分尖銳，表明算法達(dá)到收斂狀態(tài)。

5.2 算法對多目標(biāo)信號的角度分辨率

假設(shè)兩個空間方位鄰近信號分別以到達(dá)角參數(shù)θ1=-10°、θ2=-5°入射到上述傳感器陣列，陣列模型誤差百分比ρ=10%。圖2給出了采用II-MMSE算法、L1-min算法[15]和常規(guī)波束形成(CBF)算法對目標(biāo)信號到達(dá)角估計(jì)的歸一化稀疏功率譜圖。

圖2 三種不同算法得到的兩個鄰近信號到達(dá)角估計(jì)的歸一化功率譜Fig.2 Normalized power spectrums of the DOA estimate for two adjacent signals obtained by three different algorithms

由圖2中的仿真結(jié)果可看出，對于空間上方位鄰近的兩個目標(biāo)信號，采用II-MMSE算法和L1-min算法均具備對目標(biāo)信號的高方位分辨能力，而常規(guī)波束形成算法不具備這一能力。另外發(fā)現(xiàn)，L1-min算法的稀疏功率譜具有較多的小幅度偽峰，而II-MMSE算法幾乎沒有偽峰和旁瓣影響。實(shí)驗(yàn)中還發(fā)現(xiàn)，當(dāng)陣列模型誤差百分比ρ=0時，L1-min算法的偽峰也幾乎被抑制。上述結(jié)果表明，IIMMSE算法對陣列模型誤差具有一定的魯棒性，而L1-min算法不具備這一特性。

5.3 信噪比對到達(dá)角估計(jì)精度的影響

考慮兩個目標(biāo)信號分別以到達(dá)角參數(shù)θ1=-10°、θ2=60°入射到上述傳感器基陣，陣列模型誤差百分比ρ=10。對前述3種算法分別進(jìn)行300次蒙特卡洛(Monte Carlo)仿真實(shí)驗(yàn)，得到如圖3所示的目標(biāo)信號到達(dá)角估計(jì)的均方根誤差(Root Mean Square Error, RMSE)隨信噪比的變化關(guān)系。

圖3 三種算法到達(dá)角估計(jì)的RMSE隨信噪比的變化曲線Fig.3 Variation curves of the RMSE of DOA estimation with SNR by the three algorithms

從圖3中的仿真結(jié)果可發(fā)現(xiàn)，在陣列模型存在誤差的情況下，II-MMSE算法對目標(biāo)信號到達(dá)角的統(tǒng)計(jì)估計(jì)精度明顯高于另兩種算法。在低信噪比時這一優(yōu)勢更加明顯。

5.4 計(jì)算復(fù)雜度分析和計(jì)算時間比較

由式(26)可知，本文提出的II-MMSE算法的計(jì)算復(fù)雜度為O(N2M)。而L1-min算法通過文獻(xiàn)[11]中的分析可知，其計(jì)算復(fù)雜度為O(N3)。另外，由式(8)可知，CBF算法的計(jì)算復(fù)雜度為O(NM)?？紤]到在實(shí)際應(yīng)用場景中，N?M，因此在這3種估計(jì)算法中，L1-min算法的計(jì)算復(fù)雜度最高，其次是II-MMSE算法，CBF算法的計(jì)算復(fù)雜度最低。

表1中列出了上述3種估計(jì)算法在不同信噪比下的平均計(jì)算時間。實(shí)驗(yàn)中，取兩個相干信號分別以DOA參數(shù)θ1=-10°、θ2=60°入射到上述均勻線陣，陣列模型誤差百分比ρ=10%，Monte Carlo仿真實(shí)驗(yàn)次數(shù)為50。仿真環(huán)境：采用Matlab R2009a平臺，Intel Core 2處理器，2G內(nèi)存。通過比較這3種算法的計(jì)算時間可以看出，II-MMSE算法的計(jì)算效率遠(yuǎn)高于L1-min算法，但這兩種算法的計(jì)算效率都明顯低于CBF算法。但I(xiàn)I-MMSE算法和L1-min算法在整體性能上優(yōu)于CBF算法。

表1 不同信噪比下三種算法的平均計(jì)算時間Table 1 Average computation times of the three algorithms under different SNRs

6 結(jié) 論

本文通過稀疏重構(gòu)得到傳感器陣列輸出數(shù)據(jù)的稀疏表示模型，提出了一種II-MMSE算法，理論分析了II-MMSE算法的迭代收斂性和對陣列模型誤差的魯棒性，評估了該算法的計(jì)算復(fù)雜度。理論分析和仿真結(jié)果都表明II-MMSE算法保持了稀疏重構(gòu)類到達(dá)角估計(jì)方法在低信噪比、相干背景、小樣本、陣列未校準(zhǔn)等諸多非理想因素下的整體魯棒性優(yōu)勢，而且計(jì)算效率更高，無需選取正則化參數(shù)，具有潛在的工程推廣價值。