穩健的稀疏信號單快拍波達方向估計

虞 飛，宋 俊，余 赟，蘇 冰

(海軍研究院，北京 100071)

0 引 言

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

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

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

1 陣列輸出模型的稀疏重構

考慮K個遠場窄帶平面波入射到由M個無指向性單元構成的傳感器陣列，且K個信號與傳感器陣列位于同一平面內。將目標信號可能到達角范圍空間Θ進行N次均勻網格劃分，可將傳統的陣列輸出模型轉化為稀疏重構模型[11]：

2 基于MMSE準則迭代實現的單快拍DOA估計算法

根據最小均方誤差(MMSE)準則和式(1)中的稀疏重構模型，取目標函數f(t)：

式中：W(t)為M×N維復權矩陣，γ(t)顯然是經過M×N型濾波器組濾波后的期望輸出信號矢量。根據Wiener濾波器理論相關結論，容易求得式(2)中的目標函數在MMSE意義上的最優權矩陣為

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

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

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

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

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

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

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

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

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

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

式(11)、(12)和(13)構成了MMSE迭代算法第i次迭代的實現過程。當某一次迭代估計值滿足ε為某一預先指定的較小正數時，算法停止迭代，或者當算法達到預設的最大迭代次數時也停止迭代。此時，由稀疏功率對角陣的估值構成的列向量diag[(t)]關于到達角網格的稀疏功率譜中，第K個譜峰對應的角度即為目標信號到達角估計值

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

(1) 初始化

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

(2) 迭代過程

(3) 算法結果

在列向量diag[(t)]關于到達角網格的稀疏功率譜中，第k個譜峰對應的角度值即為目標信號到達角估計值

3 算法收斂性分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中：Z=diag[z1，…，zM]-IM。在式(24)的推導過程中，還利用到了Hadamard乘積的運算性質A⊙B=B⊙A和(A⊙B)H=AH⊙BH[13]。根據式(24)的結果，則權矩陣的更新變為

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

為了增強算法的自適應能力，可對式(25)進行以下修正[14]：

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

5 仿真實驗

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

實驗中，第l個陣元對應的模型誤差可表示為

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

5.1 算法收斂性實驗

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

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

從仿真結果可發現，盡管存在陣列模型誤差，II-MMSE算法仍能在真實目標方向角上形成尖銳譜峰，表明算法對目標信號具有良好的到達角估計精度和多目標分辨能力。當算法迭代到一定次數時(一般不超過15次)，功率譜圖的偽峰及旁瓣逐漸被抑制，且目標方向對應的譜峰十分尖銳，表明算法達到收斂狀態。

5.2 算法對多目標信號的角度分辨率

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

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

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

5.3 信噪比對到達角估計精度的影響

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

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

從圖3中的仿真結果可發現，在陣列模型存在誤差的情況下，II-MMSE算法對目標信號到達角的統計估計精度明顯高于另兩種算法。在低信噪比時這一優勢更加明顯。

5.4 計算復雜度分析和計算時間比較

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

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

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

6 結 論

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