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

潘曦， 佟穎， 王華陽

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

0 引言

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

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

1 基于球麥克風陣列的聲場模型

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

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

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

(1)

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

(2)

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

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

(3)

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

(4)

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

(5)

記作

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

(6)

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

(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)無關的加性噪聲向量。

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

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

(8)

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

(9)

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

(10)

Pnm(cosθ)表示締合勒讓德函數；bn(kr)稱為模態強度，且有

(11)

(12)

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

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

(13)

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

(14)

B(kr)為遠場模態強度的對角矩陣，

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)式即為聲場球諧函數展開的矩陣表達式，給(17)式左右兩邊同時乘以YH(ΩmL)×V，其中V∈CL×L表示權重因子的對角矩陣，

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

(18)

式中：αml=4π/l，其中l=1,2,…,L表示整個采樣點集。

對聲場進行球諧分解，可得球傅里葉系數矩陣為

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

(19)

式中：nnm(k)=YH(ΩmL)×V×n(k)為加性噪聲球傅里葉變換系數矩陣。注意到(19)式中球傅里葉系數向量的導向向量為B(kr)YH(ΩmL)，即頻率相關的分量已從角度相關的分量中解耦出來。這時通過給(19)式左乘B-1(kr)即可將頻率相關的分量去除，獲得只包含角度相關的分量，即球傅里葉變換成分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算法

對于寬帶聲源，其頻率覆蓋范圍寬。波數k=2πf/c的范圍由聲信號的帶寬決定，現將聲信號分解成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的行線性獨立。

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

PB=PHPA.

(25)

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

(26)

QHa(Ωs)=0，

(27)

于是DOA估計函數為

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

(28)

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

Q0=Q(QHQ)-1/2，

(29)

則DOA估計的偽譜為

(30)

通過掃描入射聲源方向Ω, 偽譜中最大的峰值對應的Ω即為聲源的入射方向Ω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)

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

CSM=(GHG)-1GHH.

(37)

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

3 仿真及性能分析

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

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

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

(38)

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

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

(39)

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

(40)

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

3.2 傳播算法DOA估計性能

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

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

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

(41)

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

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

4 實驗及測試結果分析

4.1 球陣列設計

為了驗證本文所提基于球形陣列寬帶源快速定向算法的有效性，本節設計了表面均勻分布8個全指向麥克風的球形陣列結構[14]，球陣列半徑為0.07 m，8個麥克風被安裝在剛性球內接正方體的8個頂點上，剛性球陣列實物圖如圖6所示，8個麥克風所在方位如表1所示。

表1 球陣列表面8個麥克風所在方位

4.2 定位誤差分析

將聲源固定在不同的方位進行聲源定向實驗，所用聲源信號為錄制的坦克行走聲音，該聲信號的某段時域波形如圖7所示，同時圖8給出了對應的聲源信號幅頻圖。圖8中結果顯示所用聲源頻率范圍主要集中在200～1 500 Hz之間。

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

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

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

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

4.3 時間開銷

測試了使用該單球陣列進行SH-SFT-OPM算法、SH-CSM-OPM算法和SH-MUSIC算法3種DOA估計算法花費的時間。表3給出了經過10次實驗的平均花費時間。

表3 10次實驗的平均花費時間

由表3可以看出，隨著掃描間隔的變小，3種算法DOA估計花費的時間快速地增大，其中SH-SFT-OPM算法花費的時間更少，這是因為該算法不需要計算CSM，也不需要噪聲白化。而SH-CSM-OPM算法和SH-MUSIC算法花費的時間非常相近。

5 結論

本文提出了基于SH-OPM的兩種聲源定向算法：SH-SFT-OPM算法和SH-CSM-OPM算法?？偨Y了球陣列半徑與可處理的聲源頻率范圍之間的關系，解決了寬帶聲源頻率適應性問題。通過比較分析兩種算法和SH-MUSIC算法在不同信噪比下的DOA估計性能，進行聲源定向實驗并測試了3種算法的DOA估計用時，得到了如下結論：

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

2)所提算法能夠準確識別聲源波達方向。其中SH-SFT-OPM算法估計精度較差，SH-CSM-OPM算法和SH-MUSIC算法估計結果相近。

3)隨著掃描間隔的變小，3種算法DOA估計花費的時間快速增大。其中SH-SFT-OPM算法花費的時間更少，這是因為該算法不需要計算CSM以及噪聲白化，而SH-CSM-OPM 算法和SH-MUSIC算法花費的時間相近。