基于協方差矩陣擬合的寬帶無網格波達角估計

楊益新，張亞豪，楊 龍

（1.西北工業大學航海學院，陜西西安 710072；2.陜西省水下信息技術重點實驗室，陜西西安 710072）

0 引言

寬帶波達角（Direction of Arrival,DOA）估計是聲吶系統陣列信號處理中一個非常重要的研究內容，目前比較流行的兩類算法是非相干子空間法（Incoherent Signal-subspace Method,ISM）[1]和相干子空間法（Coherent Signal-subspace Method,CSM）[2]。ISM 將陣列接收信號通過離散傅里葉變換（Discrete Fourier Transform,DFT）劃分至各個子帶，在各子帶內分別進行窄帶DOA 估計；CSM 通過聚焦矩陣將各子帶的信號轉換至聚焦頻率上，DOA 估計問題從而轉換至在聚焦頻率上的窄帶DOA 估計。與ISM相比，CSM 在低信噪比（Signal to Noise Ratio,SNR）情況下性能更加穩健。

CSM 的核心是聚焦矩陣的設計，聚焦矩陣的性能直接影響DOA 估計的精度。文獻[2-4]給出了聚焦矩陣不同的構造方法以降低聚焦誤差。這些算法需預先獲得DOA 信息作為聚焦角度來設計聚焦矩陣。因此，CSM 一般采用迭代的方式，以上一次迭代中的DOA 估計值作為下一次迭代的聚焦角度。然而，當聚焦角度與真實DOA 之間存在偏差時，將造成較大的聚焦誤差。為提升聚焦矩陣對聚焦角度選取的穩健性，陳洪光等[5]將陣列流形中的角度變量和頻率分離，僅對頻率部分進行聚焦，避免了聚焦角度的選取，然而該方法產生的聚焦誤差較大；Sellone[6]提出一種穩健聚焦矩陣，該算法可使聚焦角度周圍扇區（即聚焦區域）的聚焦誤差達到最小，因此無需得知DOA 的精確信息，僅需DOA 大致分布的先驗信息即可，從而增強了對聚焦角度估計偏差的穩健性。

完成各子帶的協方差矩陣聚焦后，現存方法一般采用子空間類算法如多重信號分類（Multiple Signal Classification,MUSIC）算法[7]和根 MUSIC（root MUSIC,rootMUSIC）[8]完成DOA 估計。然而，子空間類算法對SNR 要求較高。稀疏重構類算法[9-10]是近十幾年內發展起來的高分辨DOA 估計算法，相較于子空間類算法，這類算法可用于SNR更低的情況[9]，因此受到廣泛的關注。大部分稀疏重構類算法在DOA 估計時將空間劃分為離散的網格，一般假設信號分布在該網格上。然而，當信號偏離網格時，DOA 估計精度會有所下降。為避免該問題，近年來針對直線陣陣形發展了許多無網格類稀疏重構算法，可分為原子范數類算法[11]和協方差矩陣擬合類算法[12]。其中，原子范數類算法利用原子范數作為稀疏懲罰項，可看作為l1范數類算法向連續空間的推廣。然而，這類算法需要選取合適的正則參數控制數據擬合誤差和解的稀疏性之間的平衡；協方差矩陣擬合類算法，如稀疏與參數方法（Sparse and Parameter Approach,SPA）[12]，在協方差矩陣擬合準則下利用采樣協方差矩陣估計真實協方差矩陣，對所估計的協方差矩陣進行范德蒙分解（Vandermonde Decomposition,VD）[13]以獲取陣列流形的估計值，并從中得到DOA 估計值。相較于原子范數類算法，協方差矩陣擬合類算法無需選取正則參數，因此更易用于實際的信號處理中。

本文將SPA 算法推廣至寬帶DOA 估計中，提出一種基于 CSM 的改進 SPA（Coherent Signal-subspace based Modified SPA,C-MSPA）算法，以在低SNR 情況下獲得高DOA 估計精度。算法首先采用穩健聚焦矩陣[6]將各子帶的信號投影至聚焦頻率上，而在設計聚焦矩陣時則利用了信號空間分布的先驗信息。為提升DOA 估計對聚焦矩陣輸出的穩健性，本文基于頻率選擇的VD（Frequency-selective VD,FS-VD）原理[14]對協方差矩陣擬合準則進行改進，使得在新準則下重構的協方差矩陣中包含的DOA 信息嚴格限制在聚焦區域上，從而在DOA 估計時也有效利用了DOA 分布的先驗信息。通過仿真和湖上實測數據分析可以看出，與基于ISM 的SPA（Incoherent Signal-subspace based SPA,I-SPA）算法和其他基于CSM 的無網格類算法相比，C-MSPA算法實現了高空間方位分辨能力且提高了DOA 估計精度。

1 信號模型建立

假設K個非相關信號從θ=[θ1θ2…θK]T入射至M元均勻線列陣（Uniform Linear array,ULA），信號的中心頻率為fc，帶寬為H。將陣列接收信號劃分為N段，在每一段內進行L點的DFT，從而將時域信號轉換至頻域。在第l個子帶上，指向角度θ的陣列流形為

式中：fl和λl分別表示第l個子帶的中心頻率和對應的波長；d為陣元間距。第l個子帶上陣列接收信號可建模為

式中：xl(n) ∈?M×1，sl(n) ∈?K×1以及el(n) ∈?M×1分別表示第n段陣列接收信號、目標信號以及環境噪聲在第l個子帶上的DFT 系數；Al(θ)=[a(f l,θ1) …a(f l,θK)]∈?M×K表示陣列流形矩陣；N為頻域快拍數。

假設噪聲為均勻高斯白噪聲，且與信號不相關，則第l個子帶上陣列接收信號的協方差矩陣Rl∈?M×M為

式中：plk和σl分別為第l個子帶上第k個信號功率和噪聲功率；IM表示維度為M×M的單位矩陣；diag(x)表示以向量x為對角線元素的對角矩陣；E[?]表示期望算子。在實際信號處理中，真實的協方差矩陣一般難以獲取。因此，采用采樣協方差矩陣來近似。采樣協方差矩陣的計算公式為

2 算法原理

2.1 穩健聚焦矩陣設計

聚焦矩陣的設計原理是將頻率f l（l=1,…,L）上的信號子空間轉變為聚焦頻率f0上的信號子空間，則第l個子帶上聚焦矩陣的設計問題表示為[3]

2.2 修正稀疏與參數方法

SPA 算法的原理是在協方差矩陣擬合準則下基于采樣協方差矩陣重構真實協方差矩陣。聚焦后的協方差矩陣見式（8）。第i次迭代中聚焦后的采樣協方差矩陣表示為

基于此，SPA 算法表示為[12]

從A0(θ(i))中可很容易獲得第i次迭代中的DOA 估計值θ(i)。

3 仿真分析

下面進行仿真實驗。仿真中考慮陣元數為5 的ULA，陣元間距為0.3 m，兩個中心頻率為2.5 kHz、帶寬為1 kHz 的寬帶遠場目標信號入射至陣列。目標1 的DOA 為80.3°，目標2 的DOA 為95.2°。陣列的采樣率為104 kHz。將接收信號均勻劃分為10段，即N=10，每一段包含10 400 個快拍。每一段內進行10 400 點的DFT，將時域信號轉換至頻域，即頻率分辨率為10 Hz。分別采用C-MSPA 算法、基于 CSM 的 SPA 算法（Coherent Signal-subspace based SPA,C-SPA）、基于 CSM 的 rootMUSIC 算法（Coherent Signal-subspace based rootMUSIC,C-rootMUSIC）以及I-SPA 算法進行DOA 估計。為了保證聚焦矩陣對聚焦角度的穩健性，各基于CSM 的算法中聚焦矩陣均由式（12）計算。迭代初始的聚焦角度由CBF 算法給出，聚焦區域的初始寬度B設定為 CBF 的 3 dB 波束寬度。聚焦頻率設定為2.5 kHz，迭代終止條件設定為ε=0.01。

分別從DOA 估計精度以及分辨概率兩方面對各算法性能進行比較。其中，DOA 估計精度由均方根誤差（Root Mean Square Error,RMSE）來衡量，計算公式為

則稱兩個目標被成功分辨[17]，式中的θ1和θ2分別表示目標1 和目標2 的真實DOA。

圖1 和圖2 分別給出各算法的RMSE 和分辨概率隨SNR 的變化曲線。由圖1、2 可以看出，C-MSPA 算法在SNR 小于?10 dB 時的DOA 估計精度和分辨概率都優于其他算法。主要原因是C-MSPA 算法在進行DOA 估計時利用了信號分布的先驗信息，增強了其對噪聲的穩健性。相較于傳統的高分辨算法，稀疏重構類算法對SNR 的要求較低，C-MSPA 和C-SPA 的DOA 估計精度在SNR小于0 dB 時優于C-rootMUSIC，同時分辨概率在 SNR 小于?7.5 dB 時也高于C-rootMUSIC。I-SPA算法在各個子帶內獨立進行DOA 估計，并將各子帶的DOA 估計結果進行平均。若某一子帶的估計結果偏差較大，對該算法的估計精度會有很大影響。因此，與基于CSM 的算法相比，I-SPA 算法對SNR 的要求較高。當SNR 大于0 dB 時，I-SPA 算法的性能才慢慢接近其他算法。

圖1 各算法DOA 估計的RMSE 隨SNR 的變化 Fig.1 RMSE of each method versus SNR

圖2 各算法的分辨概率隨SNR 的變化 Fig.2 Resolution probability of each method versus SNR

為了對比算法的空間方位分辨能力，固定目標1，將目標2 從85.2°，以2°為步長移動至99.2°，即兩個目標之間的角度間隔從5°變化至19°。SNR設定為?5 dB，其他仿真條件不變。

圖3、4 分別給出了RMSE 和分辨概率隨角度間隔的變化曲線。由圖3、4 中可知，C-MSPA 算法在目標角度間隔小于9°時的RMSE 和分辨概率均優于其他算法，說明當DOA 估計時考慮信號分布的先驗信息有效提升了算法的空間方位分辨能力。由于稀疏重構類算法利用信號在空間分布的稀疏性，獲得了更高的空間方位分辨能力，因此C-MSPA 和C-SPA 算法在目標小角度間隔情況下的RMSE 和分辨概率均優于C-rootMUSIC。由圖1、2 可知，I-SPA 算法對SNR 的要求較高。由于在該仿真條件下設置的SNR 較低，導致I-SPA 算法在所有考慮的角度間隔情況下RMSE 均大于其他算法，同時其分辨概率也小于1。由此也可說明將各子帶相干處理相較于非相干處理可有效提升算法性能。

圖3 各算法的RMSE 隨角度間隔的變化 Fig.3 RMSE of each method versus angle separation

圖4 各算法的分辨概率隨角度間隔的變化 Fig.4 Resolution probability of each method versus angle separation

4 實驗分析

為了驗證C-MSPA 算法的性能，對實驗數據進行處理和分析。實驗數據來源于2016 年在新安江水庫試驗站進行的湖上實驗，湖深約為60 m。實驗設備布放如圖5 所示。

圖5 實驗設備布放圖 Fig.5 The layout of experimental equipment

實驗中，接收陣為一條自制的5 陣元ULA，陣元間距為0.3 m，水平放置在水下。實驗采用兩個聲源，其中聲源1 為低頻寬帶聲源UW350，聲源2為溢流環。根據幾何關系可以算出，聲源1 的方位角大致為75°，聲源2 的方位角大致為57°。兩個聲源同時發射中心頻率為2.5 kHz、帶寬為1 kHz 的寬帶信號，持續時間為17 s。接收陣列以104 kHz 的采樣率采集信號。將接收到的信號劃分為17 組，每一組持續時間為1 s。將每組信號劃分為10 段，即N=10，在每一段中進行10 400 點的DFT，即頻率分辨率為10 Hz。

采用第3 節仿真中的算法對每一組數據進行DOA 估計，其余的參數設置與仿真中一致，DOA估計結果如圖6 所示。考慮到CBF 算法的穩健性，圖6 同時給出了CBF 的估計結果作為參考。圖6中背景色即為CBF 的估計結果，三角符號則為各無網格算法的估計結果。

從圖6 中可以看出，由于CBF 算法的分辨能力較低，無法分辨出兩個目標，但從結果中可以看出，在區間[50，80]°存在目標。C-MSPA 算法在DOA 估計時利用了信號大致分布的先驗信息，該算法可以很好地估計出兩個目標。C-SPA 算法在DOA 估計時則沒有考慮信號的分布情況，其在某些時刻的DOA 估計值存在較大偏差。與C-MSPA和C-SPA 算法將各子帶進行相干處理不同，I-SPA算法是在各個頻點上進行非相干DOA 估計。由于部分子帶內的估計結果存在較大偏差，導致該算法無法成功分辨出兩個目標。相較于稀疏重構類算法，C-rootMUSIC 算法對SNR 要求較高，同時空間方位分辨能力弱于稀疏重構類算法，該算法也無法很好地估計出兩個目標信號。通過對實驗數據的處理，驗證了所提出的算法在寬帶DOA 估計中的優異性能。

圖6 各算法的DOA 估計結果 Fig.6 DOA estimation results of each method

5 結論

本文將協方差矩陣擬合類無網格DOA 估計算法推廣至寬帶情況。算法首先利用聚焦矩陣將各子帶的數據轉移至聚焦頻率上，以進行相干DOA 估計。不同于已有協方差矩陣擬合類算法，本文對協方差矩陣擬合準則進行改進，以利用DOA 分布的先驗信息，使在該準則下重構的協方差矩陣中所包含的DOA 信息嚴格限制在聚焦區域內。仿真和實驗數據處理結果表明，與基于ISM 的無網格算法相比，本文所提出的算法通過聚焦矩陣對各個子帶進行相干處理；相較于其他基于CSM 的無網格類算法，所提出的算法在DOA 估計時利用了DOA 分布的先驗信息，實現了較高的空間方位分辨能力且提高了DOA 估計的精度。