基于改進(jìn)MUSIC 算法的寬帶信號DOA 估計

黎昕婷，鐘舜聰，鐘劍鋒

（福州大學(xué)機(jī)械工程及自動化學(xué)院，福州 350108）

0 概述

基于傳感器陣列的波達(dá)方向（Direction of Arrival，DOA）估計是聲納、雷達(dá)、通信、語音處理等領(lǐng)域的研究熱點(diǎn)，然而現(xiàn)代化技術(shù)發(fā)展的需求使DOA 估計不再局限于處理傳統(tǒng)窄帶信號。在工程實(shí)踐中，大部分信號是非平穩(wěn)或譜時變的。在窄帶信號源假設(shè)條件下，由于源信號的導(dǎo)向矢量具有頻率不一致性，因此傳統(tǒng)子空間算法，如多重信號分類算法（Multiple Signal Classification，MUSIC）、旋轉(zhuǎn)不變子空間算法（ESPRIT）等，并不適用于工程實(shí)踐［1-3］。

近年來，研究人員對寬帶信號源條件下的DOA估計進(jìn)行了大量研究［4-6］。寬帶信號DOA 估計算法原理上［7］包含非相干信號子空間方法（Incoherent Signals-subspace Method，ISM）［8］與相干信號子空間方 法（Coherent Signals-subspace Method，CSM）［9］。非相干信號子空間方法基于頻率分解原理，將寬帶源信號分解為一系列窄帶信號后進(jìn)行DOA 估計。相干信號子空間方法基于頻率聚焦原理，其中最為典型的是雙邊相關(guān)變換［10］（Two-sided Correlation Transformation，TCT）方法，利用聚焦矩陣將分解得到的子信號變換到參考頻率點(diǎn)上，進(jìn)而使用窄帶的處理方法實(shí)現(xiàn)DOA 估計。但是該方法需要CSM 預(yù)估信號源角度，而聚焦矩陣對預(yù)估角度的依賴導(dǎo)致最終DOA 估計結(jié)果產(chǎn)生偏差。文獻(xiàn)［11］提出一種聚焦的FTOPS 算法，利用參考頻點(diǎn)的信號子空間與陣列方向矢量投影矩陣間的正交性對DOA 進(jìn)行估計。文獻(xiàn)［12］對平滑自相關(guān)矩陣進(jìn)行特征分解，再根據(jù)特征向量空間之間的過渡性構(gòu)建聚焦矩陣，從而實(shí)現(xiàn)完美聚焦的目的。但是該算法僅針對環(huán)境噪聲為高斯噪聲，并且需要預(yù)設(shè)參考頻率的情況。文獻(xiàn)［13］基于壓縮感知理論，利用陣列協(xié)方差矩陣稀疏迭代估計的方法實(shí)現(xiàn)寬帶信號DOA 估計，但是該方法在信源數(shù)目預(yù)估出現(xiàn)錯誤時，其空間譜結(jié)果易產(chǎn)生偽峰。文獻(xiàn)［14］對信號子空間聚焦法進(jìn)行改進(jìn)，采用奇異值分解重構(gòu)協(xié)方差矩陣，通過Root-正交傳播算子實(shí)現(xiàn)DOA 估計，改進(jìn)方法雖然降低了運(yùn)算量，但是仍需要預(yù)估參考頻點(diǎn)子頻帶。文獻(xiàn)［15］提出一種頻域時延補(bǔ)償方法，該方法無需對角度進(jìn)行預(yù)估，但是運(yùn)算復(fù)雜度高，難以滿足信號實(shí)時處理的要求。

針對寬帶非平穩(wěn)信號，研究人員嘗試從時頻域角度進(jìn)行研究，根據(jù)寬帶信號的時頻信息獲得更準(zhǔn)確的DOA 估計結(jié)果［16-17］?；诳烧{(diào)窗函數(shù)的S 變換和小波變換時頻分析方法能進(jìn)行多分辨率分析。文獻(xiàn)［18］構(gòu)建一種時頻域的陣列數(shù)據(jù)模型，根據(jù)信號的時頻信息來提高非平穩(wěn)信號的DOA 估計性能。文獻(xiàn)［19］利用小波包變換對信號進(jìn)行分解，再使用MUSIC 算法對每個子帶進(jìn)行空間譜估計。文獻(xiàn)［20］將S 變換應(yīng)用于MUSIC 算法，并對跳頻及交叉chirp 陣列信號進(jìn)行分析，然而該方法仍存在需要預(yù)知信號源個數(shù)的問題。文獻(xiàn)［21］提出基于小波變換的多重信號分類改進(jìn)算法，根據(jù)信號的時頻域特征來提高算法的分辨率。

本文提出一種基于改進(jìn)MUSIC 算法的寬帶信號DOA 估計。通過對接收信號進(jìn)行S 變換，獲得多分辨的時頻譜矩陣，同時構(gòu)建時頻域陣列信號模型，根據(jù)頻率段不同時刻的功率譜矩陣呈聯(lián)合對角化結(jié)構(gòu)的特點(diǎn)，設(shè)計一種新的空間譜，從而實(shí)現(xiàn)寬帶信號的DOA 估計。

1 基于S 變換的時頻陣列數(shù)據(jù)模型

假設(shè)M個陣元等間隔線性排布，陣元間距為d，已知P個寬帶信號源sp(t)從不同方向入射，入射角度分別為()θ1，θ2，…，θP，不考慮傳播過程中的信號幅值衰減，第m個陣元的接收信號如式（1）所示：

其 中：τmp=(m-1)dsinθp/v為第p個信號到達(dá)第m個陣元產(chǎn)生的時延；v為信號的傳播速度；nm(t)為第m個陣元的加性噪聲。

使用S 變換對信號xm(t)進(jìn)行處理，變換結(jié)果如式（2）所示：

其中：τ為時間因子；fi(i=1，2，…，I)為頻率分量。將上式寫成傅里葉頻譜形式進(jìn)行推導(dǎo)，如式（3）所示：

其 中：xfi，p(τ-τmp)為信號分量；nfi，m(τ)為噪聲分量。假設(shè)S 變換的頻寬和中心頻率滿足窄帶信號條件，即變換得到的信號分量為窄帶信號，如式（4）所示：

在向量形式下，陣列接收信號的時頻域模型表示如式（5）所示：

2 ST_MUSIC 算法

2.1 算法原理

根據(jù)經(jīng)典MUSIC 算法原理［22］，ST(τ，fi)的協(xié)方差矩陣如式（7）所示：

其中：RX(τ，fi)為(τ)的協(xié)方差矩陣，由于加性噪聲與源信號不相關(guān)，且自身互不相關(guān)，則變換后的噪聲方差可用表示；σi為矩陣奇異值。

由于實(shí)際接收數(shù)據(jù)長度有限，因此數(shù)據(jù)協(xié)方差矩陣取其最大似然估計，如式（8）所示：

對于第n(n=1，2，…，P)個信號源，矢量bn(f)的定義需要滿足式（9）：

不考慮噪聲項，根據(jù)功率譜矩陣的聯(lián)合對角化結(jié)構(gòu)性質(zhì)［23］，建立如下等式：

當(dāng)RY=RY(τk，f)，k=1，2，…，K時，代價函 數(shù)［23］的簡化結(jié)果如式（12）～式（14）所示：

2.2 算法步驟

本文提出的ST_MUSIC 算法主要分為以下5 個步驟：1）根據(jù)信號的有效頻段設(shè)置S 變換的頻率分量fi(i=1，2，…，I)；2）利用式（5）對接收信號進(jìn)行S 變換得ST(τ，fi)，根據(jù)式（8）計算得到協(xié)方差矩 陣3）根據(jù)式（6）構(gòu)建時頻域?qū)蛳蛄縣(θ，fi)；4）設(shè)置搜索范圍和搜索步進(jìn)Δθ，根據(jù)式（15）計算空間譜估計P(θ)；5）通過譜搜索獲得P(θ)所有譜峰所在的位置，從而得到DOA 的最大估計值。

2.3 算法分析

根據(jù)上述算法原理，ST_MUSIC 算法滿足如下條件：

1）本文算法要求符合所有寬帶信號的窄帶分量互不相關(guān)條件；

2）算法可以直接擴(kuò)展到二維源定位的情況，此時，算法的偏轉(zhuǎn)角表現(xiàn)為方位角與俯仰角的結(jié)合，平面范圍的譜搜索轉(zhuǎn)變?yōu)榭臻g譜搜索，峰值點(diǎn)坐標(biāo)為DOA 估計結(jié)果。

2.4 算法復(fù)雜度分析

本文對ST_MUSIC 算法和后續(xù)仿真中使用的雙邊相關(guān)變換方法（TCT）、基于壓縮感知理論的算法（CS_TCT）［13］及基于小波變換的MUSIC 算 法［21］（CWT_MUSIC）的計算復(fù)雜度進(jìn)行分析。本文設(shè)空間譜估計的觀測范圍的搜索點(diǎn)數(shù)為N，信號數(shù)據(jù)長為L，一個M×M維的數(shù)據(jù)協(xié)方差矩陣做特征分解，其計算量為O(M3)。文獻(xiàn)［13］已計算TCT 算法的計算復(fù)雜度為O(M3+2N3M3+N2M3)。在k個多頻點(diǎn)采用 CS_TCT 算法，其運(yùn)算復(fù)雜度為O(kN3M3+M2)。相 比TCT、CS_TCT 算 法，由 于ST_MUSIC 與CWT_MUSIC 算法中包含的I個尺度參數(shù)帶來較大的計算量，使算法復(fù)雜度顯著提升，CWT_MUSIC 算法的總計算量約為O(2IMLlbL+IM3+N)［21］。相 比CWT_MUSIC 算 法，ST_MUSIC算法在譜估計運(yùn)算中增加了O(IM)的運(yùn)算量，由于其不需要做特征值分解，總計算量仍小于CWT_MUSIC 算法，為O(2IMLlbL+N+IM)。

3 仿真結(jié)果分析

3.1 算法有效性分析

在多聲源場景下，本文對ST_MUSIC 算法進(jìn)行仿真實(shí)驗(yàn)，并與TCT、CS_TCT 以及CWT_MUSIC 這3 種算法進(jìn)行對比。本文采用16 元均勻線性陣列，構(gòu)建信噪比為0 dB、頻率范圍為165～300 Hz 的兩不相干線性調(diào)頻信號，設(shè)置信號入射角度分別為-20°和20°，采樣頻率為4 kHz，信號總長度為1 024 個數(shù)據(jù)點(diǎn)。當(dāng)陣元數(shù)為16 時，CS_TCT、TCT、CWT_MUSIC、ST_MUSIC 算法的空間譜圖如圖1 所示。從圖1 可以看出，4 種算法均可以較準(zhǔn)確地估計信號角度，但是估計效果存在差異，當(dāng)信源數(shù)估計不準(zhǔn)確時，CS_TCT 算法的偽峰抑制效果受限，而偽峰抑制效果最優(yōu)的TCT 算法在精度上較其他3 種算法略差，CWT_MUSIC 算法與本文算法具有較優(yōu)的分辨率和精度。

圖1 不同算法的空間譜圖對比（陣元數(shù)為16）Fig.1 Spatial spectrogram comparison among different algorithms（the number of array elements is 16）

為進(jìn)一步驗(yàn)證算法的有效性，本文將陣元數(shù)減少至12，其他仿真條件不變，進(jìn)行二次仿真實(shí)驗(yàn)。不同算法的空間譜圖對比如圖2 所示。從圖2 可以看出，TCT 算法在陣元數(shù)減少后出現(xiàn)多個明顯偽峰，結(jié)果出現(xiàn)錯亂，其他3 種算法對陣元數(shù)敏感度低，結(jié)果更穩(wěn)定。CS_TCT 算法估計結(jié)果準(zhǔn)確度為-20.5°和20.4°。雖然陣元數(shù)的減少不影響CS_TCT 算法最終估計結(jié)果，但是在算法仿真的零度位置出現(xiàn)較高偽峰。CWT_MUSIC 算法估計結(jié)果為-20.8°和20.9°。本文算法估計結(jié)果為-20.6°和20.6°，驗(yàn)證了本文算法的有效性。

圖2 不同算法的空間譜圖對比（陣元數(shù)為12）Fig.2 Spatial spectrogram comparison among different algorithms（the number of array elements is 12）

為進(jìn)一步探究陣元數(shù)對算法有效性產(chǎn)生的影響，在同等仿真條件下，ST_MUSIC 算法在不同陣元數(shù)下的空間譜圖如圖3 所示。從圖3 可以看出，陣元數(shù)的增加使主峰更尖銳，在提高估計精度的同時也會帶來更多的旁瓣，但其對估計結(jié)果沒有影響，而陣元數(shù)過少降低估計結(jié)果的精度，當(dāng)陣元數(shù)為4 時，DOA 估計結(jié)果偏差3°和4°，相比陣元數(shù)16，陣元估計結(jié)果的誤差保持在±0.2°。

圖3 在不同陣元數(shù)下ST_MUSIC 算法的空間譜圖Fig.3 Spatial spectrogram of ST_MUSIC algorithm under different numbers of array elements

3.2 分辨性能對比

本文采用16元均勻線性陣列，信噪比設(shè)為-10 dB，分別在信號入射角度相距120°（-45°和75°）、90°（-45°和45°）、60°（-45°和15°）、30°（-20°和10°）、10°（-5°和5°）條件下，進(jìn)行50 次隨機(jī)重復(fù)實(shí)驗(yàn)并計算4 種算法的平均分辨率。分辨率參數(shù)采用ρ=[P(θ1)+P(θ2)]/2-min{P(θ1)，P(θ1+Δθ)，…，P(θ2)}［15］，其中θ1、θ2為信號入射角，Δθ為搜索步進(jìn)，分辨率單位為dB。

在仿真實(shí)驗(yàn)中，當(dāng)信噪比降低為-10 dB 時，TCT算法的仿真結(jié)果不穩(wěn)定，并出現(xiàn)較為嚴(yán)重的偽峰，因此，本節(jié)僅對其他3 種算法進(jìn)行分辨率對比，如表1所示。從表1 可以看出，分辨率隨兩信號的入射角距增大而增高，基于壓縮感知的CS_TCT 算法相較于TCT 算法增大了角度分辨率，且分辨率結(jié)果較穩(wěn)定。本文算法在60°和120°處出現(xiàn)較高的分辨值，但是這3 種算法分辨率差距不大。

表1 不同算法的分辨率對比Table 1 Resolution comparison among different algorithms

3.3 誤差性能對比

本文仿真條件與3.1 節(jié)設(shè)置一致。本文仿真信號的信噪比范圍設(shè)置為-15～10 dB，分別通過50 次重復(fù)隨機(jī)實(shí)驗(yàn)對比4 種算法DOA 估計結(jié)果的均方根誤差，如圖4 所示。在信噪比為-15 dB 與-10 dB 時，TCT 算法估計的均方根誤差大于1°，分別為3.33°、1.60°，因此圖4 中未顯示其結(jié)果，CS_TCT 算法的DOA 估計結(jié)果在高信噪比條件下表現(xiàn)更佳。此外，在實(shí)驗(yàn)過程中，當(dāng)信噪比為-15 dB 與-10 dB 時，TCT算法的估計成功率低于50%，而CWT_MUSIC 算法與ST_MUSIC 算法在整個信噪比范圍內(nèi)估計成功率始終保持90%以上。

圖4 不同算法的均方根誤差對比Fig.4 Root mean square error comparison among different algorithms

3.4 計算量對比

本節(jié)仿真條件的設(shè)置與3.1 節(jié)保持一致。本文分別設(shè)置2 種仿真情形，分別為2 個信號（-20°和20°）和3 個信號（-60°、-20°和20°），通過50 次隨機(jī)重復(fù)實(shí)驗(yàn)對比算法的運(yùn)算時間。不同算法的平均運(yùn)算時間如圖5 所示，CS_TCT 算法的運(yùn)算時間最短，與復(fù)雜度分析結(jié)果對應(yīng)，基于頻率聚焦的TCT 算法明顯比基于頻率分解的CWT_MUSIC 算法與ST_MUSIC 算法運(yùn)算時間短，但是這4 種算法的運(yùn)算時間均在正常可接受范圍內(nèi)，ST_MUSIC 算法略優(yōu)于CWT_MUSIC 算法。

圖5 不同算法的平均運(yùn)算時間Fig.5 Average computing time of different algorithms

3.5 二維定位仿真實(shí)驗(yàn)

本文對所提算法進(jìn)行二維聲源定位仿真實(shí)驗(yàn)，其他仿真條件不變，在信噪比為0 dB 條件下，采用真實(shí)語音信號進(jìn)行實(shí)驗(yàn)，在不考慮環(huán)境混響的情況下，設(shè)置2 個語音聲源的方位角和俯仰角，分別為40°和40°、40°和-20°，定位結(jié)果如圖6 所示，從中可以看出，該算法在二維定位中能夠得到準(zhǔn)確的DOA 估計結(jié)果。

圖6 ST_MUSIC 算法的二維DOA 估計譜圖Fig.6 Two-dimensional DOA estimation spectrogram of ST_MUSIC algorithm

4 結(jié)束語

針對被動探測系統(tǒng)中的寬帶信號DOA 估計問題，本文提出一種基于S 變換且無需預(yù)估信源數(shù)的多重信號分類改進(jìn)算法。根據(jù)S 變換的多分辨率特性，通過構(gòu)建時頻陣列信號模型，提高多源DOA 估計的空間分辨率，利用譜搜索實(shí)現(xiàn)DOA 估計，實(shí)現(xiàn)多源寬帶信號的聲源定位。二維語音定位仿真實(shí)驗(yàn)驗(yàn)證了該算法的有效性。仿真結(jié)果表明，該算法具有較優(yōu)的分辨率性能和估計性能。后續(xù)將從S 變換參數(shù)的角度優(yōu)化本文算法，同時通過增強(qiáng)信號分量、弱化噪聲分量，降低算法運(yùn)算復(fù)雜度并提升其在復(fù)雜噪聲情況下的估計性能。