球諧域傳播算子快速聲定向算法

潘曦， 佟穎， 王華陽

(1.北京理工大學(xué) 機(jī)電學(xué)院， 北京 100081； 2.上海無線電設(shè)備研究所， 上海 200090)

0 引言

聲源定向是一種應(yīng)用于目標(biāo)檢測的被動定向方法。時間波達(dá)(TDOA)方法算法簡單,是普遍使用的聲源定向方法，但TDOA方法在給定的時間內(nèi)只能處理單聲源[1]定向問題，且定向精度和分辨率不高。多重信號分類(MUSIC)算法[2-4]和基于旋轉(zhuǎn)不變技術(shù)的信號參數(shù)估計(jì)(ESPRIT)算法[5-6]是真正意義上的高分辨波達(dá)方向(DOA)方法，但不能解決寬帶聲源問題[7]?；谇蛑C波理論的球形麥克風(fēng)陣列與子空間類DOA估計(jì)相結(jié)合的方法，即球諧域MUSIC (SH-MUSIC)算法能夠解決寬帶聲源問題，并具有高定向精度[8]，主要原理是將接收信號轉(zhuǎn)換到球諧域，并將球諧域的觀察空間分解成噪聲子空間和信號子空間，利用子空間的特性進(jìn)行DOA估計(jì)。然而，該方法進(jìn)行特征值分解、獲取特征向量時需要很大計(jì)算量。與MUSIC算法相似，正交傳播算子 (OPM)算法也是一種基于子空間的方法，但其不需要對傳感器接收信號的交叉譜矩陣(CSM)進(jìn)行特征值分解，并且傳播算子是一個只依賴于導(dǎo)向向量的線性算子。因此與MUSIC算法相比，OPM算法可以顯著降低計(jì)算復(fù)雜度[9]。

本文提出了基于球諧傳播算子(SH-OPM)的兩種聲源定向算法，包括使用球傅里葉變換成分獲得傳播算子(SH-SFT-OPM)的算法和使用球傅里葉變換成分的CSM獲得傳播算子(SH-CSM-OPM)的算法。實(shí)驗(yàn)結(jié)果證明，與球諧域MUSIC (SH-MUSIC)算法相比，所提基于SH-OPM的算法均具有較好的聲源定向性能和較低的計(jì)算復(fù)雜度。同時，為了實(shí)現(xiàn)對不同頻帶聲源的DOA估計(jì)，本文分析了球陣列半徑范圍從0.07～0.10 m對信號頻率范圍從0.2～2.2 kHz的估計(jì)影響問題，總結(jié)了球陣列半徑與可處理的聲源頻率范圍之間的關(guān)系。

1 基于球麥克風(fēng)陣列的聲場模型

球麥克風(fēng)陣列由球面上均勻分布的L個全指向聲壓麥克風(fēng)組成，球半徑為r，第l個麥克風(fēng)所在方向?yàn)棣竚l=(θml,φml)，在空間直角坐標(biāo)系Oxyz中第l個麥克風(fēng)的位置矢量表示為rml=(xml,yml,zml)T=(rsinθmlcosφml,rsinθmlsinφml,rcosθml)T，其中φ為方位角，θ為俯仰角，l=1, 2, …,L. 本文中所有方位角φ均在Omxmym面沿xm軸逆時針測得，所有俯仰角θ均沿zm軸向下測得。

假設(shè)有S個波數(shù)為k的寬帶源平面波入射到陣列，其中k=2πf/c，f表示聲源頻率,c表示聲速，且第s個平面波的波數(shù)為ks、第s個平面波的DOA為Ωs=(θs,φs)，s=1, 2, …,S. 那么第s個平面波的波矢量可被表示為ks=(kssinθscosφs,kssinθssinφs,kscosθs).

球陣列上第l個麥克風(fēng)接收的聲壓信號的時間- 空間域模型可以寫為

(1)

式中：ss(t)表示第s個聲源入射信號隨時間t變化的幅度；nml(t)表示麥克風(fēng)接收信號中的加性噪聲；τml(θs,φs)表示第s個聲源入射信號到達(dá)第l個麥克風(fēng)與到達(dá)球陣列球心點(diǎn)處的時間差，

(2)

Δλml為第l個麥克風(fēng)所在位置與球陣列球心點(diǎn)處在入射方向(θs,φs)上的波程差。

對(1)式進(jìn)行傅里葉變換，得到頻率- 空間域模型：

(3)

將(2)式中的τml(θs,φs)代入(3)式，得

(4)

將(4)式寫為矩陣形式：

(5)

記作

P(k)=A(k)s(k)+n(k)，

(6)

式中：A(k)∈CL×S表示導(dǎo)向矩陣, 可寫為

(7)

s(k)∈CS×1且s(k)=[s1(k),s2(k),…，sS(k)]T表示信號波形向量；n(k)∈CL×1且n(k)=[n1(k),n2(k),…，nL(k)]T表示和信號向量s(k)無關(guān)的加性噪聲向量。

與窄帶聲源信號不同，對于寬帶聲源信號的導(dǎo)向矩陣A(k)，其頻率f表示整個信號頻帶，這里自變量波數(shù)k=2πf/c不是固定值，而是整個頻帶對應(yīng)的不同波數(shù)。通過對接收信號進(jìn)行離散采樣獲得數(shù)字信號，然后選取X個快拍進(jìn)行X點(diǎn)離散傅里葉變換，就會得到X個離散的頻率點(diǎn)，相當(dāng)于將寬帶聲源信號的整個頻帶分成X個子頻帶，每個子頻帶對應(yīng)的頻率為fx(x=1, 2, …,X)，相應(yīng)的波數(shù)為kx=2πfx/c. 那么(6)式離散化為

P(kx)=A(kx)s(kx)+n(kx).

(8)

對于單位幅度的平面波入射聲場球諧函數(shù)展開式可寫為

(9)

式中：k為波矢量；r為麥克風(fēng)的位置矢量；Ynm(θ,φ)表示n階m維度的球諧函數(shù)，且

(10)

Pnm(cosθ)表示締合勒讓德函數(shù)；bn(kr)稱為模態(tài)強(qiáng)度，且有

(11)

(12)

將(12)式代入(7)式中，導(dǎo)向矩陣A(k)可被分解為

A(k)=Y(ΩmL)B(kr)YH(ΩS).

(13)

式中：Y(ΩmL)∈CL×(N+1)2為球諧函數(shù)矩陣，

(14)

B(kr)為遠(yuǎn)場模態(tài)強(qiáng)度的對角矩陣，

B(kr)=diag(b0(kr),b1(kr),b1(kr),…,bN(kr))；

(15)

Y(ΩS)∈CS×(N+1)2為入射聲源的球諧波矩陣，

(16)

將(13)式代入(8)式，得

P(k)=Y(ΩmL)B(kr)YH(ΩS)s(k)+n(k).

(17)

(17)式即為聲場球諧函數(shù)展開的矩陣表達(dá)式，給(17)式左右兩邊同時乘以YH(ΩmL)×V，其中V∈CL×L表示權(quán)重因子的對角矩陣，

V=diag(αm1,αm2,…,αml,…,αmL)，

(18)

式中：αml=4π/l，其中l(wèi)=1,2,…,L表示整個采樣點(diǎn)集。

對聲場進(jìn)行球諧分解，可得球傅里葉系數(shù)矩陣為

Pnm(k)=B(kr)YH(ΩmL)s(k)+nnm(k)，

(19)

式中：nnm(k)=YH(ΩmL)×V×n(k)為加性噪聲球傅里葉變換系數(shù)矩陣。注意到(19)式中球傅里葉系數(shù)向量的導(dǎo)向向量為B(kr)YH(ΩmL)，即頻率相關(guān)的分量已從角度相關(guān)的分量中解耦出來。這時通過給(19)式左乘B-1(kr)即可將頻率相關(guān)的分量去除，獲得只包含角度相關(guān)的分量，即球傅里葉變換成分Pnm∈C(N+1)2×1寫為

Pnm(k)=YH(ΩS)s(kx)+N(kx)，

(20)

式中：N(kx)∈C(N+1)2×1，且N(kx)=B-1(kr)×YH(ΩmL)×V×n(kx)。定義矩陣F(kx)∈C(N+1)2×L，令

F(kx)=B-1(kr)×YH(ΩmL)×V，

(21)

則N(kx)可被重寫為

N(kx)=F(kx)×n(kx).

(22)

2 球諧域的傳播算子定向算法

2.1 SH-SFT-OPM算法

對于寬帶聲源，其頻率覆蓋范圍寬。波數(shù)k=2πf/c的范圍由聲信號的帶寬決定，現(xiàn)將聲信號分解成X個子頻帶，則球傅里葉變換成分Pnm(k)寬帶信號的球傅里葉變換成分Pnm∈C(N+1)2×X可被寫為

Pnm=[Pnm(k1),Pnm(k2),…,Pnm(kX)].

(23)

傳播算子的定義基于對球傅里葉變換成分Pnm獲得

Pnm=[PA;PB],

(24)

式中：PA∈CS×X表示由矩陣Pnm的前S行組成的子矩陣；PB∈CU×X表示由矩陣Pnm的后U行組成的子矩陣，U=(N+1)2-S.

假定所有的入射信號均不相干，則矩陣Pnm滿秩保證了子矩陣PA的行線性獨(dú)立。

在無噪聲的情況下，存在著矩陣P使得(25)式成立：

PB=PHPA.

(25)

當(dāng)存在噪聲時，盡管(24)式仍然成立，(25)式關(guān)系已不再滿足。這時傳播算子SFT估計(jì)值可以通過最小二乘法獲得，可表示為

(26)

QHa(Ωs)=0，

(27)

于是DOA估計(jì)函數(shù)為

FPM(Ωs)=aH(Ωs)QQHa(Ωs).

(28)

為了引入噪聲空間的投影算子，可以用其正交化矩陣取代矩陣Q，即

Q0=Q(QHQ)-1/2，

(29)

則DOA估計(jì)的偽譜為

(30)

通過掃描入射聲源方向Ω, 偽譜中最大的峰值對應(yīng)的Ω即為聲源的入射方向Ωs.

2.2 SH-CSM-OPM算法

X個子頻帶的CSMRnm∈C(N+1)2×(N+1)2可被寫為

(31)

(32)

(33)

(34)

(35)

在無噪聲的情況下，存在著矩陣P使得(36)式成立：

H=GP.

(36)

當(dāng)存在噪聲時，盡管(24)式和(35)式的分塊仍然成立，但(25)式和(36)式的關(guān)系已不再滿足。這時傳播算子CSM的估計(jì)值可以通過最小二乘法獲得，可表示為

CSM=(GHG)-1GHH.

(37)

獲取傳播算子后，構(gòu)建矩陣Q以及DOA估計(jì)都與SH-STF-OPM算法相同。

3 仿真及性能分析

在基于球形麥克風(fēng)陣列的DOA估計(jì)中，DOA的估計(jì)性能主要取決于球傅里葉變換成分Pnm的準(zhǔn)確性。本節(jié)通過定義球傅里葉變換成分的信噪比(SNR-SFT)對球麥克風(fēng)陣列的DOA估計(jì)性能進(jìn)行分析。

3.1 球陣列半徑與可處理的聲源頻率范圍

為了考慮不同來源的誤差，通過隨機(jī)過程仿真球形麥克風(fēng)陣列的性能。給定L個麥克風(fēng)以一個給定的安裝誤差均勻地分布在球形陣列的表面。用球陣列上L個麥克風(fēng)接收到的聲壓信號計(jì)算S個平面波的來波方向[12]，且所有的接收聲壓值都添加了隨機(jī)噪聲。使用球形陣列的SNR-SFT可寫為

(38)

式中：‖·‖F(xiàn)表示Frobenius范數(shù)；nm(k)表示球傅里葉變換成分的估計(jì)值。Pnm(k)和nm(k)可分別寫為

Pnm(k)=YH(Ωs)s(k)，

(39)

nm(k)=F(k)A(k)s(k)+F(k)n(k).

(40)

結(jié)合圖1和圖2可以得出：在低頻處球傅里葉變換成分的估計(jì)值nm并不真正等于其理論值Pnm，這種現(xiàn)象控制了球麥克風(fēng)陣列可處理頻率范圍的下邊界；在高頻處由于空間混疊的存在，限制了球麥克風(fēng)陣列的性能[13]。對半徑和陣列可處理的頻率范圍作了相應(yīng)的分析，結(jié)果如圖3所示。從圖3中可以看出，隨著球陣列半徑的增加，陣列DOA估計(jì)可處理的聲源信號頻率范圍先增加、后下降。

3.2 傳播算法DOA估計(jì)性能

本節(jié)所有的仿真均使用同一個遠(yuǎn)場正弦線性掃頻信號，該信號頻率范圍為0.2～2.2 kHz，并以8 kHz/s頻率變化率線性掃頻。聲信號的入射方向(θ,φ)為(100°, 100°)。麥克風(fēng)接收信號的信噪比設(shè)置為0 dB. 快拍數(shù)為512個，采樣率為10.24 kHz.

使用半徑為0.07 m剛性單球陣列，進(jìn)行仿真實(shí)驗(yàn)，圖4給出了3種DOA估計(jì)算法的空間譜。通過比較可以得出：圖4(a)中使用SH-SFT-OPM 算法的DOA估計(jì)準(zhǔn)確性最差；圖4(b)和圖4(c)中分別使用SH-CSM-OPM算法和SH-MUSIC算法的DOA估計(jì)結(jié)果非常接近，和圖4(a)中使用SH-SFT-OPM 算法得到的空間譜相比，這兩種算法的空間譜譜峰均很尖銳。

通過比較3種算法DOA估計(jì)結(jié)果的均方根誤差(RMSE)分析這3種算法在不同信噪比下的DOA估計(jì)性能。RMSE由 (41)式計(jì)算：

(41)

式中：S和W分別表示聲源個數(shù)和獨(dú)立的蒙特卡洛實(shí)驗(yàn)次數(shù)；(s,w,s,w)表示第w次蒙特卡洛實(shí)驗(yàn)中以(θs,φs)方向入射到陣列第s個聲源的DOA估計(jì)值[14]。俯仰角和方位角的掃描間隔設(shè)為(1°, 1°)，麥克風(fēng)接收信號的信噪比為-10～20 dB. 入射聲源信號的所有參數(shù)即3.2節(jié)開頭所設(shè)。使用半徑為0.07 m的剛性單球陣列對SH-SFT-OPM算法、SH-CSM-OPM算法和SH-MUSIC算法3種DOA估計(jì)算法進(jìn)行估計(jì)誤差分析，共進(jìn)行500次獨(dú)立蒙特卡洛實(shí)驗(yàn)，得到RMSE隨信噪比(SNR)變化的結(jié)果如圖5所示。

從圖5可知：本文所提SH-SFT-OPM算法和SH-CSM-OPM 算法在信噪比為0～20 dB時具有與SH-MUSIC算法相似的性能；在信噪比低于0 dB時所提算法的RMSE比SH-MUSIC算法小，特別是SH-SFT-OPM算法不需要計(jì)算CSM，也不需要進(jìn)行白化，因此該算法具有最低的計(jì)算復(fù)雜度，在對實(shí)時性要求較高的環(huán)境中具有實(shí)際應(yīng)用價值。

4 實(shí)驗(yàn)及測試結(jié)果分析

4.1 球陣列設(shè)計(jì)

為了驗(yàn)證本文所提基于球形陣列寬帶源快速定向算法的有效性，本節(jié)設(shè)計(jì)了表面均勻分布8個全指向麥克風(fēng)的球形陣列結(jié)構(gòu)[14]，球陣列半徑為0.07 m，8個麥克風(fēng)被安裝在剛性球內(nèi)接正方體的8個頂點(diǎn)上，剛性球陣列實(shí)物圖如圖6所示，8個麥克風(fēng)所在方位如表1所示。

表1 球陣列表面8個麥克風(fēng)所在方位

4.2 定位誤差分析

將聲源固定在不同的方位進(jìn)行聲源定向?qū)嶒?yàn)，所用聲源信號為錄制的坦克行走聲音，該聲信號的某段時域波形如圖7所示，同時圖8給出了對應(yīng)的聲源信號幅頻圖。圖8中結(jié)果顯示所用聲源頻率范圍主要集中在200～1 500 Hz之間。

圖9為實(shí)際聲源方向測量實(shí)驗(yàn)場景，用藍(lán)牙音箱分別在距離球陣列中心2 m處8個不同的空間位置播放聲音，圖10給出了8個聲源空間位置的示意圖，圖中灰色小球表示球陣列，赤道圓環(huán)上的4個黑色小點(diǎn)表示實(shí)驗(yàn)位置1到實(shí)驗(yàn)位置4中俯仰角固定為90°時方位角分別為0°、90°、180°、270°的位置，4個藍(lán)色小點(diǎn)表示實(shí)驗(yàn)位置5到實(shí)驗(yàn)位置8中固定方位角為180°時，俯仰角分別為30°、60°、120°、150°的位置。

每個聲源位置都進(jìn)行了100次獨(dú)立實(shí)驗(yàn)，并使用3種不同算法進(jìn)行計(jì)算，所得結(jié)果的平均值如表2所示，使用SH-CSM-OPM算法進(jìn)行第3次測量結(jié)果的空間譜三維圖和俯視圖分別如圖11(a)和圖11(b)所示。

表2 8個不同空間位置測得的聲源方向平均值

從表2中可以看出，所提算法能夠準(zhǔn)確識別聲源DOA. 其中SH-SFT-OPM算法估計(jì)精度較差，最大估計(jì)誤差為(4°, 6°)；而SH-CSM-OPM算法和SH-MUSIC算法估計(jì)結(jié)果相近，最大估計(jì)誤差為(3°, 5°)。綜上所述，實(shí)驗(yàn)結(jié)果證明了所提算法的有效性。

4.3 時間開銷

測試了使用該單球陣列進(jìn)行SH-SFT-OPM算法、SH-CSM-OPM算法和SH-MUSIC算法3種DOA估計(jì)算法花費(fèi)的時間。表3給出了經(jīng)過10次實(shí)驗(yàn)的平均花費(fèi)時間。

表3 10次實(shí)驗(yàn)的平均花費(fèi)時間

由表3可以看出，隨著掃描間隔的變小，3種算法DOA估計(jì)花費(fèi)的時間快速地增大，其中SH-SFT-OPM算法花費(fèi)的時間更少，這是因?yàn)樵撍惴ú恍枰?jì)算CSM，也不需要噪聲白化。而SH-CSM-OPM算法和SH-MUSIC算法花費(fèi)的時間非常相近。

5 結(jié)論

本文提出了基于SH-OPM的兩種聲源定向算法：SH-SFT-OPM算法和SH-CSM-OPM算法。總結(jié)了球陣列半徑與可處理的聲源頻率范圍之間的關(guān)系，解決了寬帶聲源頻率適應(yīng)性問題。通過比較分析兩種算法和SH-MUSIC算法在不同信噪比下的DOA估計(jì)性能，進(jìn)行聲源定向?qū)嶒?yàn)并測試了3種算法的DOA估計(jì)用時，得到了如下結(jié)論：

1)所提算法在信噪比為0～20 dB時具有和SH-MUSIC算法相似的性能；在信噪比低于0 dB時所提算法的RMSE比SH-MUSIC算法小，特別地，SH-SFT-OPM算法具有最低的計(jì)算復(fù)雜度，在對實(shí)時性要求較高的環(huán)境中具有實(shí)際應(yīng)用價值。

2)所提算法能夠準(zhǔn)確識別聲源波達(dá)方向。其中SH-SFT-OPM算法估計(jì)精度較差，SH-CSM-OPM算法和SH-MUSIC算法估計(jì)結(jié)果相近。

3)隨著掃描間隔的變小，3種算法DOA估計(jì)花費(fèi)的時間快速增大。其中SH-SFT-OPM算法花費(fèi)的時間更少，這是因?yàn)樵撍惴ú恍枰?jì)算CSM以及噪聲白化，而SH-CSM-OPM 算法和SH-MUSIC算法花費(fèi)的時間相近。