使用同心多環陣提升聲源定位魯棒性?

張國昌 吳 鳴 韓欣宇 武帥兵 殷 蘭 楊 軍

(1中國科學院大學 北京 100049)

(2中國科學院噪聲與振動重點實驗室聲學研究所 北京 100190)

(3北京信息科技大學自動化學院 北京 100192)

0 引言

近年來，使用麥克風陣列的聲源方位估計方法引起了學者的廣泛關注，其在機器人、電話會議、視頻監控、助聽器等[1-4]領域有著諸多的應用。通常聲源波達方向(Direction of arrival,DOA)估計方法分為間接法和直接法兩類[5-7]。間接法，該方法需要首先計算聲源到各陣元對的相對時間延遲(Time difference of arrival,TDOA)而后使用幾何方法確定聲源方位[8]；直接法，其計算一系列候選方位上的損失函數，并選擇最有可能的方位，其根據使用的技術不同分為可控響應功率(Steered response power,SRP)法[9]、最大似然方法[10]和子空間類方法[11]。

其中TDOA方法雖然計算復雜度較低，但是由于其只利用了各通道信號互相關函數最大值所對應的時延信息，使得其極易受混響的影響。

由廣義互相關函數(Generalized cross correlation,GCC)計算SRP的DOA估計方法被廣泛的應用分析，其中基于相位變換(Phase transform,PHAT)廣義互相關函數的SRP算法(Steered response power-phase transform,SRP-PHAT)[12-15]最為流行。但是這種方法的空間分辨率較差，另外還需要進行逐個頻點的空間譜掃描，計算復雜度高。

最大似然方法以其良好的一致性(Consistency)、極佳的統計特性(Statistical efficiency)和處理相干信號的能力被廣泛的研究[10]，但需要進行高維空間搜索來求取最優解，計算復雜度高，目前很難實際應用。

另外一種基于子空間的方法以其超分辨特性也受到了廣泛的關注，這類方法在多聲源或存在混響的情況下具有優于TDOA和SRP方法的性能[11]。針對寬帶信號，兩種基于窄帶分解的子空間方法得到廣泛的研究。非相干信號子空間方法(Incoherent signal subspace method,ISSM)[8]融合各個子帶的DOA估計結果；相干子空間方法(Coherent signal subspace method,CSSM)[16-19]先將各窄帶信號聚焦至某一特定頻點，而后只在聚焦頻點處使用窄帶DOA估計。兩種方法各有優缺點：前者由于相干早期反射聲的存在使得空間相關矩陣(Spatial correlation matrix,SCM)的信號子空間缺秩，這導致估計性能的急劇下降；后者利用頻率聚焦解決了前者的缺秩問題，但需要對聲源方位進行初步估計以獲得聚焦矩陣，該方法的性能對初始DOA估計精度非常敏感[18]。

為了避免初步DOA估計和聚焦矩陣的實時求解，針對環形陣列(Circular array,CA)，環諧波(Circular harmonic,CH)分解可以解耦頻率與估計方位角之間的相關性，使得聚焦可以通過各頻點下SCM的平均來實現。本文后續的聲源定位方法均采用CH域的CSSM。

CH域的DOA估計方法以其優良的特性被廣泛的分析。Tianaroig等[20]使用環諧波波束形成(Circular harmonic beamforming,CHB)改善了延時求和波束形成器(Delay and sum beamformer,DSB)的分辨力；Teutsch等[21]提出了特征波束域的基于旋轉不變技術的信號參數估計方法(Estimating signal parameter via rotational invariance techniques,ESPRIT)，避免了空間譜掃描；Torres等[22]首次使用逐時頻點處理的方法來提升CHB定位的穩健性。但是CA的CH域展開系數存在零點，這會導致嚴重的噪聲放大問題。為了解決這一問題，以下三種方法被廣泛應用：將CA固定至剛性圓柱散射體上，使用吉洪諾夫正則化[23]和使用同心多環陣(Concentric circular array,CCA)[24]。與前兩種方法相比，多環陣無需增加剛性散射體也不存在正則化引入的波束失真。

為了解決CH展開零點處的噪聲放大問題，本文由CH域的理想波束出發，推導出了CCA頻域補償濾波器的最小模求解方法。這種方法設計的補償濾波器具有最優的噪聲增益。另外，最小模CCA環諧波域的波束響應也不存在理想CA貝塞爾函數零點處的失真和正則化CA的低頻波束展寬問題。為了對比CA和最小模CCA在CH域的聲源DOA估計性能，還設計了一套麥克風陣列系統。仿真和實驗均表明，與使用同陣元數、同孔徑的CA相比，最小模CCA可以有效地降低聲源方位估計的誤差。

1 同心多環陣環諧波展開

1.1 連續同心多環陣

圖1展示了由P個連續同心單環陣組成的連續同心多環陣列，其中第p個環陣的半徑為rp，信源s入射的方位角和俯仰角分別為θi和φi。第p個連續環陣接收到的聲壓信號可在極坐標下表示為

其中，S(ω)為入射信源信號在角頻率ω下的分量，k=2πf/c=ω/c為波數，c為聲音在空氣中的傳播速度。

圖1 多環陣信號入射示意圖Fig.1 The geometry of CCA with P rings

利用貝塞爾函數展開公式

式(2)中，Jl(α)為一類l階貝塞爾函數，式(1)可以展開成

這 里ejlθ被 稱 為l階 環 諧 波CH(或 模 態)，Cl(krp,θi,φi)為由(θi,φi)方向入射信源信號的ω頻率成分在l階CH上的展開系數，

至此，第p個連續圓環陣列在各方位角下的聲壓信號被分解到無窮階的CH上，xp(krp,θ)與Cl(krp,θi,φi)為一傅里葉變換對。

1.2 連續陣列均勻采樣

實際應用時，需要對環形連續陣列進行離散采樣，這里考慮均勻采樣場景，這種離散的單環陣和同心多環陣列分別被稱為均勻環形陣列(Uniform circular array,UCA)和均勻同心環形陣列(Uniform concentric circular array,UCCA)。采樣操作引入了CH系數的誤差[23,25]：

其中，為采樣引起的第p環l階CH系數誤差，假設第p環的CH分解階數為lp，那么當第p環麥克風數目Np＞|2lp|時，式(5)中的第一項占主要成分。另外CH分解階數lp＞krpsin(φi)時的CH分解系數較小可以忽略。考慮sin(φi)取最大值1，p環的CH最高階數假定為Lp≈krp，由于CH最高階數為krp，所以當第p環麥克風數滿足Np＞2krp時較小可以被忽略。

2 多環陣環諧波域波束形成

多環陣CH域波束形成旨在將各環各階CH成分進行合適的補償，以形成具有空間篩選特性的波束。其關鍵在于補償系數的求解。下面將從理想波束出發，推導出UCCA的最優補償濾波器。

理想情況下，針對入射聲源s，指向θi方向的期望波束形成器的響應輸出可表示為

這里δ(θ-θi)可以由無窮階的傅里葉級數組合獲得，則

CHB的輸出可以表示成各環各階環諧波分量進行合適補償后的累加，

為了使得CHB的輸出可以逼近于式(7)所示的理想波束形成器的輸出，應當滿足以下條件：對任意給定環數P，存在一組或者多組可使得式(10)成立：

這里ωL至ωH為設計補償濾波器的頻率范圍。式(10)中的在環數P=1時有唯一確定解，在P＞1時有無窮個可行解。

考慮均勻離散采樣，第p環CH分解最高階數Lp，并假設L1=L2=···=Lp=L，理想波束近似為

將式(11)寫成向量形式得

此時，BCHB(ω,θ)可以表示如下：

CH域波束BCHB(ω,θ)可以看作是CH系數矢量和CH域導向矢量的夾角余弦。

2.1 單環陣補償濾波器求解及其正則化

當環數P=1時，UCCA退化為UCA，此時補償濾波器Hl(ω,φi)需要滿足

在頻率ω對應貝塞爾函數零點時，該頻點對應的Hl(ω,φi)趨于無窮大，這導致了嚴重的噪聲放大問題。為了緩解這一問題，一種基于吉洪諾夫正則化求解補償濾波器的方法被提出[23]

式(16)中，α為吉洪諾夫正則化系數，正則化的引入提升了CHB的穩健性，但犧牲了波束指向性。增大正則化系數可以獲得更加穩定的波束形成器。

2.2 多環陣最小模補償濾波器求解

多環的引入，松弛了式(10)對應的約束條件，使得補償濾波器由UCA的唯一解變為UCCA的無窮個解。可以從這無窮個解中挑選出最魯棒的解。與陣元域的加權系數的特性類似，最魯棒解對應的補償濾波器趨向于具有最小的模值[26]。最小模補償濾波器可以通過求解如下最優化問題獲得：

Hl(ω,φi)=為各環補償濾波器的向量形式，

式(8)對應的最小模解為

與UCA利用標量Hl(ω,φi)進行CH系數補償不同，UCCA使用矢量濾波器同時補償多個環陣的同階CH。為了補償某一確定(l,ω,φi)條件下的CH系數至期望值，UCCA可以利用不同半徑下的多個同階CH來完成，從而具有更大的靈活性。通過求解最小模補償濾波器，在接近貝塞爾函數零點處，UCCA會偏向于賦給CH分解系數較小的環陣以較小的補償系數，這一點與式(16)對應的UCA補償濾波器的解恰恰相反。通過合理的設置半徑，可以避免各環CH分解系數零點處于相同頻率條件下，UCCA可以徹底解決UCA中貝塞爾函數零點處的噪聲放大問題。

補償濾波器模值與CH波束形成器的穩健性是負相關的[26]。下面對UCA和UCCA在俯仰角為90°下的前兩階補償濾波器的模值進行了仿真，其中，正則化參數配置為0.00065，兩種UCA和UCCA的陣列參數如表1所示。結果如圖2所示，可知，UCA低階CH下存在模值無窮大的補償系數(圖中展示的是模值的倒數)，UCCA通過組合兩個無公共貝塞爾函數零點的UCA完全解決了這一問題。另外值得注意的是，與沒有進行正則化處理的兩UCA相比，UCCA在任一頻點下都具有更小的補償系數模值，也就是說，UCCA在任一頻率下都具有更加穩健的性能。圖3對比了正則化后UCA和UCCA的波束響應，其中，CH分解階數為3階，單環陣半徑選擇為60 mm。正則化雖然有效地提升了單環陣CHB的穩健性，但并不能有效地消除貝塞爾函數零點導致的波束圖失真。正則化UCA的波束響應在2180 Hz、3440 Hz、4630 Hz等頻點處存在明顯的畸變，另外，正則化還引起了低頻波束的展寬。而UCCA則不存在上述問題。由于CH分解階數的限制，在高頻區域，DSB方法具有最窄的主瓣。

表1 陣列參數Table1 Configurations of microphone arrays

圖2 兩種理想UCA、正則化UCA和UCCA補償濾波器模值與頻率的關系曲線Fig.2 The compensation filters’modulus values curves are shown as a function of frequency for UCAs with no regularization,UCA with regularization,and UCCA

圖3 單環DSB、CHB和同心雙環CHB的波束圖Fig.3 Beampattern for DSB,CHB using UCA with regularization and CHB using UCCA

3 DOA估計

3.1 逐時頻區域的構建CH域SCM

與傳統UCA構建SCM的方法類似[27]，UCCA的CH域SCM可以由式(19)計算得到：

3.2 使用MUSIC方法進行DOA

在每個時頻區域，DOA都會被使用多重信號子空間分類(Multiple signal classification,MUSIC)的方法估計出來。在假設每個TF區域只存在一個聲源的前提下，聲源方位角可以通過搜索空間譜最大值來獲得：

這里A(θ)為方位角θ對應的導向矢量，UN∈C(2L+1)×2L為由SCM最小的2L個特征值對應的特征矢量組成的噪聲子空間，Θ為待掃描的方位區間。

4 仿真及實驗

4.1 仿真

本小節，為了驗證UCCA具有更優的定位魯棒性，該文評估了不同信噪比和混響條件下的DOA估計的統計性能。評估算法均采用CH域的CSSM方法。聲源語料來自于TIMIT[28]語料庫。UCA半徑和陣元數分別設置為60 mm和16個，UCCA配置參數如表1所示，其他的仿真參數設置如表2所示。UCA的最高CH展開階數為7階，UCCA的CH展開階數受限于陣元數最少的環，最高為3階。在接下來的性能評估實驗中，麥克風陣列被放置于距離地面1 m高的房間水平面中心位置，單個播放聲源的揚聲器被依次放置于方位角間隔30°的6個方位，揚聲器高度同為1 m，聲源到麥克風陣列中心位置距離為2 m。以下仿真結果會統計平均所有聲源方位。虛源法[29]被用于生成帶有混響的多通道信號，多通道不相關的高斯白噪聲被添加至帶混響信號，加入噪聲后的信號被看作是陣列采樣信號。另外，需要使用語音激活檢測(Voice activity detection,VAD)模塊[30]對陣列采樣信號進行預處理，截取語音段，舍去非語音段。語音段數據每10幀(0.16 s)獲得一個平均的聲源方位估計結果。每一幀數據選取離散傅里葉變換后的88個子帶(500～6000 Hz)，考慮頻域步進間隔為5個子帶，此時每一幀數據會進行16次聲源方位估計。那么10幀一共獲得160個定位結果，對該160個結果進行中值平均進而估計出一個較為穩健的聲源方位。仿真分析中使用的軟件為Matlab 2018b，處理器為Intel Core i7-6700HQ，主頻2.6 GHz，處理器的核心數為4。

表2 仿真參數設置Table2 Parameters of simulations

首先仿真不同信噪比條件下的UCA和UCCA的DOA估計的統計性能，將房間的混響時間(T60)設置為0.5 s，信噪比依次設置為5 dB、10 dB、15 dB和20 dB。仿真結果如表3所示，其中最后一行的實時計算系數的計算公式為實時系數等于聲源定位算法的執行時長除以VAD檢測后的數據時長。

由表3可知，UCA在CH展開階數為3階或4階的時候取得了較優的平均絕對離差性能。過低或者過高的分解階數均導致了魯棒性的退化，這是因為過低的分解階數會使得過多的空間信息被舍棄，過高的分解階數又會引起較差的白噪聲增益。UCCA在階數為3階的時候具有較優的魯棒性。對比UCA和UCCA，僅僅在高信噪比(20 dB)CH展開階數為4階這一個條件下，UCA具有優于UCCA的性能。在該條件下，因為UCCA構建的CH域空間相關矩陣維度較低，所以具有更低的計算復雜度。

接下來仿真不同混響時間下的DOA估計性能，設定信噪比為15 dB，混響時間設置范圍為0.2～1.4 s，絕大部分房間的混響時間也均處在這一范圍。表4展示了不同混響時間下的聲源定位的統計性能。由表4可知，當T60＞0.2 s時，3階CH展開下的UCCA具有最優的魯棒性，當T60為0.2 s時，3階CH展開的UCCA與4階CH展開的UCA具有接近的性能。另外，與UCA相比，相同CH展開階數的UCCA的魯棒性具有顯著的優勢。這也就是說，在計算復雜度相同的情況下，UCCA性能優勢明顯。

表3 不同信噪比和陣列參數下的DOA估計平均絕對離差Table3 Mean absolute deviation performance for various SNRs and different array configurations(單位：°)

表4 不同混響時間和陣列參數下的DOA估計平均絕對離差Table4 Mean absolute deviation performance for various T60s and different array configurations(單位：°)

值得注意的是，目前UCCA各環的半徑設置只遵循以下準則：沒有同階CH展開系數的零點處在相同頻率點下。這一設計準則不是最優的，有理由認為，存在某種最優準則下的環半徑設計方法可以使得UCCA具有更優的聲源定位性能。

4.2 實驗

仿真未考慮實際使用場景下的麥克風一致性、陣元位置誤差等因素，因此設計一套UCCA用于評估其在真實環境下的定位魯棒性具有重要意義。采用樓氏電子微機電系統(MEMS)模擬麥克風設計一套參數如表1所示的UCCA麥克風陣列硬件，如圖4所示。麥克風型號為SPH1642，其幅度靈敏度誤差為±1 dB，信噪比為65 dB。使用該陣列，在一個經過裝修的房間中進行實驗。房間尺寸為5.8 m×4.6 m×2.4 m。實驗時揚聲器和麥克風陣列布置方法如圖5所示，每次實驗只有一個揚聲器發聲音。其中揚聲器和麥克風陣列的放置高度均為距離地面1.2 m，揚聲器距離麥克風陣列中心1.7 m。

圖4 含有16個麥克風的均勻同心雙環陣實物圖，半徑分別為60 mm、40 mmFig.4 Photograph of 16-element uniform concentric circular microphone array with radius of[60 mm,40 mm]

圖5 實驗配置Fig.5 Experimental setups

在實驗開始前，首先將揚聲器放置于圖5所示的0°方向測量了房間的脈沖響應，并使用脈沖反向積分方法估計了房間的T60，房間實物圖和房間脈沖響應如圖6所示，脈沖反向積分法[31]估計出T60約為0.377 s。

圖6 實驗房間實物圖及其對應的房間脈沖響應Fig.6 Photograph of the experiment room and room impulse response of this room

接下來選用與仿真相同的參數和語料數據進行了實驗。實驗時依次將揚聲器放置在圖5所示的7個位置，每個方位播放30 min的語料。麥克風陣列采樣信號首先進行VAD處理，處理后的數據每0.16 s獲得一次中值平均的定位結果。表5展示了不同CH展開階數和不同聲源方位下DOA估計的統計性能。表格中的UCA對應為UCCA的外環。

表5 不同聲源位置和陣列參數下的DOA估計平均絕對離差Table5 Mean absolute deviation performance for various source localizations and different array configurations(單位：°)

與仿真結果類似，在CH展開階數相同時，UCCA的定位魯棒性在所有實驗布放的聲源位置下均具有明顯優勢。分析CH展開階數，在大部分的聲源方位下，UCA在CH展開階數為4階時定位性能較優，但與3階CH展開的UCCA相比，仍有不小的差距。

綜合上述的仿真和實驗結果來看，UCCA在以下條件下均具有顯著優于 UCA的定位性能：(1)相同CH展開階數，上述所有的T60和SNR條件下；(2)強混響或者低信噪比時，任意的CH展開階數下。在高噪聲和弱混響的部分CH展開階數下，UCA與UCCA具有接近的DOA估計性能。

5 結論

本文闡述了一種針對多環陣列的環諧波展開技術。環諧波展開系數經補償濾波器補償后被用于聲源方位估計。為了提升CH域方位估計的穩健性，針對UCA貝塞爾函數零點處噪聲放大問題，提出了一種最小模準則設計UCCA補償濾波器的方法。針對混響環境，利用環諧波展開解耦頻率與方位角相關性的特點，CSSM被用于UCA和UCCA。仿真結果表明，合理設置孔徑的最小模UCCA可以解決UCA零點問題，同時也避免了正則化引入的低頻波束展寬。統計性能的仿真和實驗顯示，在相同麥克風數目和陣列孔徑的前提下，與UCA相比，使用UCCA可以顯著提升混響環境下的DOA穩健性。