基于協(xié)方差矩陣擬合的寬帶無網(wǎng)格波達(dá)角估計(jì)

楊益新，張亞豪，楊 龍

（1.西北工業(yè)大學(xué)航海學(xué)院，陜西西安 710072；2.陜西省水下信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，陜西西安 710072）

0 引言

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

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

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

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

1 信號模型建立

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

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

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

假設(shè)噪聲為均勻高斯白噪聲，且與信號不相關(guān)，則第l個(gè)子帶上陣列接收信號的協(xié)方差矩陣Rl∈?M×M為

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

2 算法原理

2.1 穩(wěn)健聚焦矩陣設(shè)計(jì)

聚焦矩陣的設(shè)計(jì)原理是將頻率f l（l=1,…,L）上的信號子空間轉(zhuǎn)變?yōu)榫劢诡l率f0上的信號子空間，則第l個(gè)子帶上聚焦矩陣的設(shè)計(jì)問題表示為[3]

2.2 修正稀疏與參數(shù)方法

SPA 算法的原理是在協(xié)方差矩陣擬合準(zhǔn)則下基于采樣協(xié)方差矩陣重構(gòu)真實(shí)協(xié)方差矩陣。聚焦后的協(xié)方差矩陣見式（8）。第i次迭代中聚焦后的采樣協(xié)方差矩陣表示為

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

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

3 仿真分析

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

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

則稱兩個(gè)目標(biāo)被成功分辨[17]，式中的θ1和θ2分別表示目標(biāo)1 和目標(biāo)2 的真實(shí)DOA。

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

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

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

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

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

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

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

4 實(shí)驗(yàn)分析

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

圖5 實(shí)驗(yàn)設(shè)備布放圖 Fig.5 The layout of experimental equipment

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

采用第3 節(jié)仿真中的算法對每一組數(shù)據(jù)進(jìn)行DOA 估計(jì)，其余的參數(shù)設(shè)置與仿真中一致，DOA估計(jì)結(jié)果如圖6 所示?？紤]到CBF 算法的穩(wěn)健性，圖6 同時(shí)給出了CBF 的估計(jì)結(jié)果作為參考。圖6中背景色即為CBF 的估計(jì)結(jié)果，三角符號則為各無網(wǎng)格算法的估計(jì)結(jié)果。

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

圖6 各算法的DOA 估計(jì)結(jié)果 Fig.6 DOA estimation results of each method

5 結(jié)論

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