一種基于重要抽樣的麥克風陣列時延估計

摘 要： 麥克風陣列信號處理已經廣泛引用于語音處理、通信等領域。時延的精確估計決定了麥克風陣列聲源定位系統的性能。將重要抽樣的方法應用到麥克風陣列的時延估計中，可以保證時延估計精確性的前提之下，有效地回避了網格搜索和迭代計算，使結果不依賴于初始值，并且降低了計算的復雜度。仿真實驗結果表明，重要抽樣時延估計方法在計算量降低一個數量級時，仍保持著與ML網格搜索法相近的性能，更適合麥克風陣列時延的實時估計。

關鍵詞： 麥克風陣列； 聲源定位； 時延估計； 重要抽樣

中圖分類號： TN911?34 文獻標識碼： A 文章編號： 1004?373X（2015）24?0036?04

Time delay estimation of microphone array based on important sampling

YANG Bo1， LIU Zhicheng1， LI Demei2

（1. Engineering University of CAPF， Xi’an 710086， China； 2. Inner Mongolia PAP Corps of CAPF， Hohhot 010000， China）

Abstract： Microphone array signal processing has been widely used in the fields of speech processing and communication. The performance of microphone array sound source localization system is determined by the accurate estimation of time delay. If the important sampling method is applied to the time delay estimation of microphone array， the grid search and iterative calculation can be effectively avoided in the precondition of guaranteeing the accuracy of time delay estimation， so that the results do not depend on the initial value， and the computational complexity is reduced. The simulation results show that the time delay estimation method based on the important sampling keeps the performance similar to ML grid search method while the calculated amount is reduced by an order of magnitude， and is more suitable for the real?time estimation for the time delay of the microphone array.

Keywords： microphone arrays； acoustic source localization； time delay estimation； important sample

麥克風陣列的聲源位置估計的問題。實質上就是對陣元之間接收信號時延的估計；因此，能否準確和快速的估計出時延，決定了麥克風陣列是否能夠實時的、準確的估計出聲源的位置。時延的估計在很多領域當中都是一個重要問題，例如雷達、聲納和無線通信網絡[1]。傳統方法采用匹配濾波和互相關的技術來進行時延估計。眾所周知的最大似然時延估計ML方法，在時延估計方面展示出了良好的性能，尤其在高信噪比和信號特性緩慢變化時可以方差逼近著名的克拉美羅界（CRLB）下限。盡管如此，ML方法在計算的復雜度上仍存在一定缺陷，基于網格搜索的ML估計計算復雜度會隨著參數數量增加而急劇上升。當然，采用迭代的方法成功的解決了ML的問題，但是估計結果的準確度過分依賴于初始值，如果沒有一個好的初始值，迭代的結果就會出現較大偏差[2]。采用ML方法對多個未知參數進行估計時，引入重要抽樣[3]的概念可以有效地避開網格搜索和迭代所產生的計算復雜度問題。

1 系統模型

假設麥克風陣列總共含有P個陣元，每個陣元都將接收信號傳輸到一個中心節點進行處理，麥克風陣列的接收信號可以表示為[4]：

[yt=i=1Pαixt-τi+w（t）] （1）

式中：[xt]是聲源信號；[αi]是衰減因子，即第i個麥克風的接收信號相對聲源信號的衰減系數，在[0，1]區間內取值（近場模型中，[αi]各不相同；遠場模型中，[αi]均相同，可以歸一化為1）；[τi]為第i個陣元接收信號同聲源信號之間的時延；[w（t）]是加性噪聲，為了簡化模型，假設不同陣元間噪聲相互獨立與信號不相關。

在接收端，以[Fs]的采樣頻率對信號進行采樣，采樣率[Fs=1Ts]，采樣信號可以表示為：

[ynTs=i=1PαixnTs-τi+wnTs， n=0，1，2，…，N-1] （2）

式中N為樣本總數。

將接收端樣本轉換到頻域，模型就可以用矩陣的形式表達。對式（2）進行離散傅里葉變換，可以得到：

[Yk=i=1PαiXkexp-j2πkτiN+Wk] （3）

其中：[{Y（k）}N-1k=0]，[{X（k）}N-1k=0]和[{W（k）}N-1k=0]分別是[{y（nTs）}N-1k=0]，[{x（nTs）}N-1k=0]和[{w（nTs）}N-1k=0]的離散傅里葉變換。

將第1個陣元當作參考陣元，其接收信號采樣的頻域表示為：

[Y1k=α1Xkexp-j2πkτiN+W1k] （4）

那么陣列接收信號與參考陣元接收信號之間的時延[Δτi=τi-τ1]，i=1，2，…，P，顯然[Δτ1=0]。為了突顯出待估計的參數，對式（3）進行處理得到：

[Yk=i=1PβiY1kexp-j2πkΔτiN+Wpk] （5）

式中：

[βi=αiα1， i=1，2，…，P] （6）

[Wpk=Wk-βiW1（k）exp-j2πkΔτiN] （7）

然而式（5）滿足一般線性模型（GLM），可以將其轉化為如下形式：

[Y=ΦpΔτβ+Wp] （8）

式中：矩陣[Y]是麥克風陣列所接收到的信號；矩陣[ΦpΔτ]僅取決于未知的時延[Δτ=[Δτ1，Δτ2，…，][Δτi…，Δτp]]。

[ΦpΔτ=ΦpΔτ1，ΦpΔτ2，ΦpΔτi，…，ΦpΔτp] （9）

向量[{ΦpΔτi}Pi=1]可以定義成：

[ΦpΔτi=Y10，Y11e-j2πΔτiN，Y12e-j2π2?τiN，…， Y1N-1e-j2π（N-1）?τiNT] （10）

其中：[Y1=[Y10，Y11，Y12，…，Y1N-1]T]是與參考麥克風接收脈沖相一致包含傅里葉系數的樣本，[Wp=[Wp0，Wp1，Wp2，…，WpN-1]T]是加性噪聲成分。在這里并沒有必要保證噪聲嚴格服從高斯分布，噪聲樣本從連續時間過程得到之后，根據中心極限定理，每一個傅里葉系數[W（k）]都可以看作對這些樣本的一個簡單加權，如果接收樣本數量足夠多，噪聲就可以被模型化為一個高斯隨機變量。那么，似然函數[5]可以表示如下：[ΛΔτ，β∝pY；Δτ，β= 1πMσ2MWexp -1σ2W×Y-Φp（Δτ）βH（Y-Φp（Δτ）β）] （11）

式中：[pY；Δτ，β]是[Δτ]和[β]參數化的概率密度函數；[σW]是[Wk]的方差，即噪聲功率。

傳統的ML估計方法是通過對式（11）搜索似然函數最大時所對應的[Δτ]來得到[ΔτML]。然而，[ΛΔτ，β]由[Δτ]和[β]共同決定，進行它們的聯合估計計算復雜度較大。因此，[有必要得出一個僅由未知量Δτ]決定的似然函數。為此可以考慮，將[β]表示為與[Δτ]相關的形式[β（Δτ）]。[β（Δτ）]通過由一個給定的[Δτ]來對對數似然函數[LΔτ，β=log {ΛΔτ，β}]搜索最大值所對應的[β]得到。

將式（11）中[β]用[β（Δτ）]替換，然后將與[Δτ]無關的部分略去，就得到了如下的系統壓縮似然函數：[LcΔτ=1σ2WYHΦpΔτ（ΦpΔτHΦpΔτ）-1ΦpΔτHY] （12）

2 時延的估計

理論上通過找出[LcΔτ]最大值所對應的[Δτ]即可估計出來未知時延。若直接用多維網格搜索的方法來求解這個問題，計算復雜度也會隨著未知參數數量的增加呈指數方式增加。一些用來解決多維最優化問題的可選途徑展示出了相對好的結果，但是大部分途徑涉及到了迭代運算。因此需要有一個良好的初始值，否則就無法保證收斂于全局最優解。在這個情況之下，Pincus提出了一個可以確保得到多維函數全局最優解的方法。根據[6]中提到的定理，代價函數[LcΔτ]的多維全局最優解[Δτ=[Δτ1，Δτ2，…，][Δτi…，Δτp]T]，可以表示如下：

[Δτi=limρ→∞∫J…∫JΔτiexp {ρLc（Δτ）}dΔτ∫J…∫Jexp{ρLc（Δu）}dΔu] （13）

式中J是待估計時延[Δτ]被限定的取值區間。

定義一個偽概率密度函數：

[L‘c，ρ0Δτ=exp {ρ0LcΔτ}∫J…∫Jexp{ρ0Lc（Δu）}dΔu] （14）

那么，[Δτi]的最優值為：

[Δτi=∫J…∫JΔτiL‘c，ρ0ΔτdΔτ] （15）

當[ρ0]取值趨近無窮時，[L‘c，ρ0Δτ]就變為了一個P維狄拉克（Dirac）函數。它的極大值所在位置會趨近[LcΔτ]的極大值所在位置，即使[Δτ]并不是真正的隨機變量，[L‘c，ρ0Δτ]包含了概率密度函數的所有性質。 然而式（15）需要計算多維積分，直接進行求解存在很大困難，不難想到[Δτi]是[Δτi]的估計值，是依據[L‘c，ρ0Δτ分布的Δτ]當中的第i個元素。那么，就可以用蒙特卡洛（Monte?Carlo）方法有效估計出[Δτ]的期望值[Δτ]，可以表示如下：

[Δτ=1Rk=1RΔτk] （16）

[{Δτk}Rk=1]是根據[L‘c，ρ0Δτ]生成向量[Δτ]的R次實現。這樣就可以把式（15）中的多重積分用一個簡單的樣本平均替換。此時新的問題又出現了，[L‘c，ρ0Δτ]相對于待估計參數[Δτ]是非線性的，并不適用于簡單的生成[Δτ]。這時需要構造與[L‘c，ρ0Δτ]相似并且簡單易于生成實現[Δτ]的函數。重要抽樣方法能夠有力地解決多重積分的問題，觀察如下的關系式：

[J…JfΔτL‘c，ρ0ΔτdΔτ=J…JfΔτL‘c，ρ0Δτg′（Δτ）g′（Δτ）dΔτ ] （17）

式中：[g′（Δτ）]是另一個偽概率密度函數，稱為歸一化重要性函數（IF）；[f·]是一個給定的函數。

如此可得式（17）左邊就是[fΔτL‘c，ρ0Δτg'（Δτ）]的均值，其中[Δτ]根據[g′（Δτ）]生成。因此，在構造[g′（Δτ）]的時候，既要保持與[L‘c，ρ0Δτ]的相似性，又要保證易于生成[Δτ]。根據蒙特卡洛（Monte?Carlo）方法，可以將式（17）寫成如下形式：

[∫J…∫JfΔτL‘c，ρ0Δτg′Δτg′ΔτdΔτ=1Rk=1RfΔτkL‘c，ρ0Δτkg′（Δτk）] （18）

其中[Δτk]是根據[g′Δτ]對[Δτ]的第k次實現。

由式（12）可以看出，[LcΔτ]主要由[（ΦpΔτHΦpΔτ）-1]構成，為了避免矩陣的轉置，可以用對角陣[k=1NY1k2-1Ip]來近似表示[（ΦpΔτHΦpΔτ）-1]。[Ip]是[p×p]確定矩陣。這樣IF函數[gρ1（Δτ）]可以定義為：

[gρ1Δτ=expρ1σ2Wk=1NY1k2YHΦpΔτΦpΔτHY]

[=i=1Pexp ρ1σ2Wk=1NY1k2I（Δτi）] （19）

式中：[ρ1]是常數；[I（Δτi）]是在頻域求[Δτi]值的數據周期圖，如下所示：

[IΔτi=k=1NY1k*Y（k）exp j2π（k-1）ΔτiN2] （20）

將[gρ1Δτ]歸一化處理得到歸一化重要性函數（IF）：

[gρ1′′Δτ=i=1Pexp {ρ1′I（Δτi）}∫J…∫Ji=1Pexp {ρ1′I（Δui）}dΔui] （21）

[ρ1′=ρ1σ2Wk=1NY1k2] （22）

理論上歸一化的IF函數[gρ1′′Δτ]的圖像含有P個峰，每個峰出現在不同時延所對應的位置。因此每實現一次[gρ1′′Δτ]，就會生成一組[Δτk]，[Δτk（i）]分別在每一個峰所對應的位置附近生成。然而由于噪聲影響，會出現一些次級波峰對生成時延造成影響。增大參數[ρ1′]可以減少次級波峰，但是同時也可能毀壞有用的波峰。最優的[ρ1′]是當IF函數剛好出現P個波峰的時候，所對應的最大的[ρ1′]值。

式（21）呈現出很多優勢，比如可以將每個時延在[gρ1Δτ]中的聯合貢獻分離成了每個時延的個體貢獻的乘積形式[7]。因此，通過[gρ1Δτ]對[Δτ]進行實現也就可以轉換成對[Δτ]中的每一個[Δτi]的實現，同時也就將求[Δτ]的問題轉化成求[Δτi]的問題。

根據（18）可以得到：

[Δτi=1Rk=1RΔτk（i）FΔτk] （23）

[FΔτk=L‘c，ρ0Δτkgρ1′′Δτk] （24）

根據式（23）可以估計出[Δτ=[Δτ1，Δτ2，…，Δτi…，Δτp]T]。

[FΔτk]作為一個加權系數，盡可能減小偏離[Δτi]的實現對估計[Δτi]時產生的影響。為了避免指數運算可能造成的計算數值溢出，可以將[FΔτk]用[F′Δτk]進行替換[8]：

[F′Δτk=expρ0LcΔτk-ρ1′i=1PIiΔτki- max1≤l≤R（ρ0LcΔτl-ρ1′i=1PIiΔτli）] （25）

由于待求參數[Δτi]具有上界和下界，這里用三角均值的概念來求解[Δτi]比較合適[9]，因此，可以表述成如下的形式：

[Δτi=12πL∠1Rk=1RF（Δτk）exp j2πΔτk（i）L] （26）

式中L是待估計時延可能的取值區間長度。

3 仿真和討論

這里采用計算機的仿真來檢驗上述的算法，仿真使用Matlab語言。

仿真過程中，聲源為遠場聲源，入射角設置為與法線夾角45°，采用語音信號。在室溫下，設置3陣元間距為0.1 m的均勻直線麥克風陣列，麥克風采樣率為44.1 kHz，信號幀長為4 096個采樣點。第一個陣元設為參考陣元，求陣元和參考陣元接收信號的時延。為了更客觀地對重要抽樣的時延估計方法進行評估，這里將其與ML的網格搜索方法進行比較。

圖1、圖2描繪出了聲源是語音信號時，時延估計在不同信噪比（SNR）下的均方誤差（MSE）。觀察實驗結果，重要抽樣的方法的性能可以趨近ML網格搜索法，當信噪比大于0時，可以估計出較為準確的時延，在高信噪比時，性能也處于同一數量級，但是在計算復雜度上比ML網格搜索法低一個數量級。

圖1 估計第一個時延的均方誤差（MSE）

圖2 估計第二個時延的均方誤差（MSE）

對算法的復雜度進行一個簡要的分析。若需要估計P個時延，每個時延進行N個點的搜索，對[LcΔτ]每1次搜索要進行5次矩陣乘法運算、2次矩陣共軛轉置、1次矩陣求逆。對[gρ1Δτ]每1次搜索要進行2次矩陣共軛轉置、3次矩陣乘法。那么，網格搜索法就對[LcΔτ]進行了[NP]次搜索，重要抽樣法對[gρ1Δτ]進行了[P×N×R]次搜索，R為對[Δτ]的實現生成次數。隨著P的增加，R也要相應增大才能保證時延估計得性能不下降。設N=100，P取1～4，R對應取10，50，100，250。圖3描繪出了估計時延個數不同時，兩種算法的大致計算復雜度對數值。顯然隨著估計時延個數的增加，ML網格搜索法計算復雜度急劇上升，而重要抽樣法計算復雜度緩慢上升。因此，重要抽樣法更適合于麥克風陣列時延的實時估計。

圖3 計算復雜度對比

4 結 語

將重要抽樣的思想應用到麥克風陣列的時延估計當中，是對傳統的極大似然估計ML方法的一個改進，有效避免了多維網格搜索和迭代運算所造成的計算復雜度以及初始值的問題。相比于傳統的方法，在保證時延估計準確的前提下，有效降低了計算復雜度，更適合于時延的實時估計。

參考文獻

[1] 劉建平，張一聞，劉穎.基于麥克風陣列的汽車笛語識別及笛聲定位方法[J].西安電子科技大學學報，2012，39（1）：163?167.

[2] HAHN P J， MATHEWS V J， TRAN T D. Adaptive realization of a maximum likelihood time delay estimator [C]// the 1996 IEEE International Conference on Acoustics， Speech， and Signal Processing. Atlanta， GA： IEEE， 1996： 3121?3124.

[3] 陶巍，劉建平，張一聞.基于麥克風陣列的聲源定位系統[J].計算機應用，2012，32（5）：1457?1459.

[4] NG L， BAR SHALOM Y. Multisensor multitarget time delay vector estimation [J]. IEEE Transactions on Acoustics， Speech and Signal Processing， 1986， 34（4）： 669?677.

[5] PINCUS M. A closed form solution for certain programming problems [J]. Operations Research， 1968， 16（3）： 690?694.

（上接第39頁）

[6] PALLAS M， JOURDAIN G. Active high resolution time delay estimation for large BT signals [J]. IEEE Transactions on Signal Processing， 1991， 39（4）： 781?787.

[7] SAHA S， KAY S. An exact maximum likelihood narrowband direction of arrival estimator [J]. IEEE Transactions on Signal Processing， 2008， 56（4）： 1082?1092.

[8] MCKILLIAM R G， QUINN B G. Estimating the circular mean from correlated observations [M]. [S.l.]： Daspworkshop Org， 2011.