穩(wěn)健的稀疏信號(hào)單快拍波達(dá)方向估計(jì)

虞 飛，宋 俊，余 赟，蘇 冰

(海軍研究院，北京 100071)

0 引 言

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

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

雖然稀疏重構(gòu)類(lèi)到達(dá)角估計(jì)方法在整體魯棒性方面優(yōu)于子空間類(lèi)估計(jì)技術(shù)，但是前者的主流算法在計(jì)算復(fù)雜度方面處于相對(duì)劣勢(shì)，而且還涉及額外的正則化參數(shù)選取問(wèn)題。本文通過(guò)稀疏重構(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é)方差信息，對(duì)于陣列模型誤差具有較好的魯棒性。II-MMSE算法保持了稀疏重構(gòu)類(lèi)到達(dá)角估計(jì)方法在諸多非理想因素下的整體魯棒性?xún)?yōu)勢(shì)，克服了主流稀疏重構(gòu)類(lèi)方法的高計(jì)算復(fù)雜度問(wèn)題，而且避免了正則化參數(shù)選取，進(jìn)一步推動(dòng)了稀疏重構(gòu)類(lèi)到達(dá)角估計(jì)方法的工程實(shí)用化。

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

考慮K個(gè)遠(yuǎn)場(chǎng)窄帶平面波入射到由M個(gè)無(wú)指向性單元構(gòu)成的傳感器陣列，且K個(gè)信號(hào)與傳感器陣列位于同一平面內(nèi)。將目標(biāo)信號(hà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)過(guò)M×N型濾波器組濾波后的期望輸出信號(hào)矢量。根據(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é)方差矩陣，對(duì)于功率為σ2n的空間高斯白噪聲，Rn可以簡(jiǎn)化為Rn=σ2nIM，P=E[γ(t)γH(t)]為N×N維稀疏信號(hào)協(xié)方差矩陣。假設(shè)稀疏信號(hào)矢量γ(t)內(nèi)各個(gè)分量互不相關(guān)，則有：

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

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

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

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

則稀疏功率對(duì)角陣的初始值為

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

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

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

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

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

(1) 初始化

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

(2) 迭代過(guò)程

(3) 算法結(jié)果

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

3 算法收斂性分析

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

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

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

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

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

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

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

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

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

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

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

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

其中：Δ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維零均值矢量。

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

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

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

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

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

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

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

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

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

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

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

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

圖1 II-MMSE算法迭代不同次數(shù)得到的兩個(gè)信號(hào)到達(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)方向角上形成尖銳譜峰，表明算法對(duì)目標(biāo)信號(hào)具有良好的到達(dá)角估計(jì)精度和多目標(biāo)分辨能力。當(dāng)算法迭代到一定次數(shù)時(shí)(一般不超過(guò)15次)，功率譜圖的偽峰及旁瓣逐漸被抑制，且目標(biāo)方向?qū)?yīng)的譜峰十分尖銳，表明算法達(dá)到收斂狀態(tài)。

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

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

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

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

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

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

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

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

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

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

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

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

6 結(jié) 論

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