虞 飛, 余 赟, 周利輝, 彭春光
(中國人民解放軍92578部隊, 北京 100071)
信號波達方向(direction of arrival,DOA)估計是陣列信號處理中一項極為重要的研究分支,已在聲納、雷達、通信、導航、地震探測等領域取得廣泛應用[1]。以MUSIC[2-3]和ESPRIT[4]算法為代表的傳統子空間類DOA估計算法在高信噪比(signal to noise ratio, SNR)、非強相關源和多采樣快拍數的情況下可獲得超分辨估計性能。但在實際應用中,高SNR、非強相關源和多采樣快拍數都是比較理想的環境因素,只要有一個因素不滿足,這些算法的DOA估計精度將嚴重下降甚至估計失敗。為此,學者們提出了很多改進算法,如適用于單快拍的陣列測向方法[5-7],適用于相干情形的空間平滑類算法[8-10]等。然而,這些改進算法都只適用于某一個非理想環境因素下的高性能估計,當同時存在至少兩個非理想環境因素時,這類算法依然面臨著失效問題。
近年來,基于稀疏表示框架的陣列參數估計方法由于在低SNR、有限采樣快拍數據、相干信號、信號DOA角度間隔小等非理想條件下同時具有良好的估計性能而引起了相關學者的廣泛關注,并取得了大量高質量的研究成果[11-15]。考慮到在實際應用中,通常目標信號個數遠小于傳感器陣元數,目標信號DOA相對于空間來說也是稀疏的,故可將傳統的傳感器陣列輸出模型進行稀疏化表示,得到陣列輸出數據的稀疏表示模型。于是,稀疏表示理論中的很多方法很自然地被引入到陣列信號處理問題中,而傳統的陣列參數估計問題轉化為稀疏參數求解問題。基于稀疏表示的陣列測向方法在非理想條件下整體性能優于傳統的子空間類測向方法,但分析表明,前者的計算復雜度明顯高于后者,導致算法對目標信號DOA估計的實時性變得很差。同時,這類方法在求最優解過程中,還面臨著一個甚至多個超參數(或者正則化參數)選取的問題,超參數的選擇是否合適直接關系到最終的稀疏恢復性能。因此,還需要對超參數進行精細選取,然而目前還沒有明確的方法用來指導超參數的選擇。
本文依托陣列輸出數據的稀疏表示模型,提出了一種基于加權最小二乘(weighted least squares,WLS)準則的單快拍DOA稀疏迭代自適應估計(簡稱為WLS-IAE)算法。該算法不僅保持了稀疏表示框架在非理想條件下的良好估計性能,而且避免了傳統稀疏參數求解方法的高復雜度和經驗式的超參數選取問題,僅采用單快拍即獲得高精度DOA估計,因此非常適用于快變目標信號DOA的實時跟蹤測量,具有潛在的工程實用價值。
考慮K個遠場窄帶信號入射到由M個無方向性(全向)陣元構成的均勻線性陣列(uniform linear array,ULA)中,并假設K個信號與ULA在同一平面內(在實際應用場景中,如水下聲納基陣測向,我們經常只關心某一平面內信號的入射方位,或者入射信號在該平面內的投影,因此這一假設可以得到保證)[16]。將陣元由1到M進行編號,并以陣元1作為基準或參考陣元。設參考陣元處的任一接收信號變換到基帶后有如下形式:
s(t)=a(t)ej[ω0t+φ(t)]
(1)

那么,在t時刻,整個陣列的M×1維輸出數據模型為
x(t)=As(t)+n(t)
(2)
式中,s(t)=[s1(t),s2(t),…,sK(t)]T為變換到基帶后的參考陣元處接收到的K個信號構成的列向量;n(t)=[n1(t),n2(t),…,nM(t)]T為陣列的零均值加性復高斯白噪聲向量;A=[a(θ1),a(θ2),…,a(θK)]為M×K維導向矢量矩陣,且對于ULA,導向矢量a(θk)[17]可定義為
(3)
式中,d為相鄰陣元間距;λ為信號波長;θk為第k個信號的DOA,通常定義為該信號入射方向與ULA法線方向的夾角,則有θk∈[-π/2,π/2]。

x(t)=Φγ(t)+n(t)
(4)


min‖γ(t)‖0s.t.‖x(t)-Φγ(t)‖2≤β
(5)
式中,l0-范數‖γ(t)‖0表示向量γ(t)中非零元素的個數;‖·‖2表示向量的2-范數;正則化參數β用來指定允許的噪聲水平。
直接求解優化問題式(5),必須篩選出向量γ(t)中所有可能的非零元素,由于搜索空間過于龐大,故此方法是非確定性多項式時間(non-deterministic polynomial-time, NP)難題。
考慮到l1-范數是最接近于l0-范數的凸目標函數,目前,使用最廣泛的求解方法是將式(5)的l0-范數最小化問題轉化為凸松弛的l1-范數最小化問題(簡稱l1-min算法)[18],即
min‖γ(t)‖1s.t.‖x(t)-Φγ(t)‖2≤β
(6)
式中,‖·‖1表示向量的l1-范數。該式可以很容易轉化為二階錐規劃(second-order cone programming,SOCP)問題的一般形式,從而在二階錐規劃的框架下求解。但在實際應用時,目標可能來向的角度范圍θ在網格劃分時其數量級將達到104,這將明顯降低二階錐規劃問題的求解速度,導致算法對DOA估計的實時性變得很差。

假設稀疏信號矢量γ(t)與噪聲矢量n(t)之間相互統計獨立,則由式(4)可得

ΦPΦH+Rn
(7)

P=diag[p1,p2,…,pn,…,pN]
(8)

定義干擾(即稀疏信號矢量γ(t)中除了第n個信號γn(t)之外的其他所有信號)的協方差矩陣為
(9)

(10)

為求得式(10)的最優解,將該式展開并處理得如下形式:
(11)
式中,
故關于γn(t)的加權最小二乘估計為
(12)
將式(9)代入式(12),根據矩陣求逆引理可得
(13)

(14)

(15)
相應地有
(16)

重復計算:
forn=1,2,…,N
end
直至收斂。


(17)


(18)

(19)
式中,向量em表示M×M維單位陣IM=[e1,e2,…,eM]的第m列。
綜上,WLS-IAE算法最終形式歸納如下。
初始化:
重復計算:
forn=1,2,…,N
end
直至收斂。
l1-min算法和本文提出的WLS-IAE算法都是基于稀疏表示框架的單快拍DOA估計方法,其中,l1-min算法作為依賴于超參數的稀疏估計類方法具有一定代表性,在其基礎上發展了很多改進算法。下面詳細分析WLS-IAE算法的計算量,并與l1-min算法進行比較。
令MDN表示復數乘除法運算的次數。首先分析本文算法初始化的計算量,再分析算法迭代過程的計算量。

(20)
結合WLS-IAE算法初始化過程的計算量,則WLS-IAE算法實現所需的MDN總次數約為
(21)
l1-min算法一般是將式(6)的凸松弛的約束優化問題轉化為二階錐規劃問題的一般形式,并在二階錐規劃的框架下采用內點法求解,其計算復雜度為O(N3)[18]。
考慮到在實際應用中,M?N,nitr≈10,則nitr?N,且一般有nitr?M成立,則WLS-IAE算法的計算復雜度為O(M2N)。
因此,雖然同為基于稀疏表示框架的DOA估計算法,WLS-IAE算法的計算復雜度遠低于l1-min算法,具有更好的實時性。
在以下仿真實驗中,考慮由15個陣元構成的均勻線性傳感器陣列,相鄰兩個陣元間距為信號波長的一半,即d=λ/2。假設兩個單位功率的窄帶信號分別從不同方向入射到上述傳感器陣列,實驗中取陣列對信號的單快拍采樣數據。SNR定義為
下列實驗中,如無特別說明,取SNR=10 dB。
(1)算法收斂特性
考慮兩個相干信號DOA參數分別以θ1=-10°,θ2=60°入射到上述均勻線陣。圖1給出了采用WLS-IAE算法分別在初始化、迭代1次、迭代5次、迭代10次過程中對DOA估計的歸一化稀疏功率譜圖。其中算法初始化采用的是CBF算法。兩個相干信號之間的關系表示為s2(t)=s1(t)?ejη,其中,?ejη為相關系數,反映了兩個信號之間的數值關系。實驗中取?=1,η=π/6。

圖1 WLS-IAE算法中DOA估計的歸一化稀疏功率譜
從圖1的4組仿真曲線可以看出,本文提出的WLS-IAE算法在相干信號的真實目標方向上形成了譜峰。隨著迭代次數的逐漸增加,稀疏功率譜的譜峰變得越來越尖銳,而“偽峰”及旁瓣逐漸消失,意味著該算法在迭代10次時對目標信號DOA具有較高的估計精度和分辨率。在改變其他實驗參數的情況下,算法一般不超過15次迭代即可達到收斂狀態,說明本文算法具有較快的收斂速度。
(2)算法對相干信號DOA估計的分辨率
考慮兩個相關系數為ejπ/6的鄰近相干信號參數θ1=-10°,θ2=-6°入射到上述均勻線陣,采用CBF算法、WLA-IAE算法和l1-min算法對單快拍陣列接收數據分別進行10次蒙特卡羅仿真實驗,得到3種算法對DOA估計的歸一化稀疏功率譜圖,分別對應如圖2(a)~圖2(c)所示。

圖2 角度間隔較小的兩個相干信號DOA估計的歸一化稀疏功率譜
由圖2的仿真曲線可以看出,當相干目標信號的角度間隔較小時,采用WLS-IAE算法和l1-min算法均可以在真實的目標方向形成尖銳的譜峰,實現了對目標信號DOA的高分辨率估計。相反,CBF算法在兩個真實目標方向上的譜峰合二為一,說明此時CBF算法無法分辨兩個信號的DOA。綜上表明,在單快拍情形下,基于稀疏表示的陣列測向方法對空間鄰近信號仍具有較好的分辨能力。
(3)SNR對DOA估計精度的影響
考慮兩個相關系數為ejπ/6的相干信號DOA參數分別以θ1=-10°,θ2=60°入射到上述均勻線陣,采用WLA-IAE算法、l1-min算法、CBF算法、快速迭代內插波束形成(fast iterative interpolated beamformer,FIIB)算法[5]和稀疏協方差矩陣迭代的單快拍(sparse covariance matrix iteration with a single snapshot,SCMISS)算法[15]分別進行200次蒙特卡羅仿真實驗,得到目標信號DOA估計的均方根誤差(root-mean-square error,RMSE)隨SNR的關系曲線,如圖3所示。其中,L次蒙特卡羅仿真實驗得到的DOA估計的RMSE定義為

圖3 DOA估計的RMSE隨SNR的關系曲線

從圖3的仿真曲線可以看出,隨著SNR的逐漸提高,采用上述5種算法估計的DOA的RMSE都是逐漸減小的,說明5種算法對DOA的估計精度均越來越高。其中,WLA-IAE算法對DOA的估計精度最高,尤其在低SNR時估計精度優勢更明顯,表明WLA-IAE算法在相干信號、單快拍、低SNR等非理想情形下均具有很好的估計性能。
CBF算法的估計性能因受陣列孔徑的限制,顯然不如WLS-IAE、l1-min和SCMISS超分辨算法,它們都屬于基于稀疏表示框架的陣列測向方法,在低SNR情況下都具有良好的估計性能。但是,l1-min算法的稀疏恢復性能依賴于正則化參數β的精細化選取,一旦選取不當,算法存在收斂至局部極值的風險。雖然本文在第3節給出了正則化參數β的選擇依據,但正如文獻[18]所述,這種選擇依據只適用于高SNR情況下,在低SNR時對正則化參數β只能憑經驗公式選取,而WLS-IAE算法屬于一種不依賴超參數的稀疏估計類方法,因此WLS-IAE和l1-min算法在高SNR時性能相當,而在低SNR時前者性能優于后者。
本文依托陣列輸出數據的稀疏表示模型,提出了一種WLS-IAE算法。詳細分析了WLS-IAE算法的基本性能和計算復雜度,并與依賴于超參數的經典稀疏估計類算法l1-min的計算復雜度進行了比較。仿真結果表明:WLS-IAE算法在低SNR、單快拍、信號相干、目標信號空間間隔很小等非理想情況下都具有很好的估計精度和分辨率;WLS-IAE算法的計算效率明顯高于一般的稀疏估計類方法。因此,WLS-IAE算法不僅保持了稀疏表示框架在非理想條件下的良好估計性能,而且避免了傳統稀疏參數求解方法的高復雜度和經驗公式的超參數選取問題,因此非常適用于快變目標信號DOA的實時跟蹤測量,具有潛在的工程實用價值。