基于粒子群優(yōu)化的TDOA聲源定位方法

張大桂,周志峰,張 怡,王立端

(1.上海工程技術(shù)大學(xué) 機(jī)械與汽車工程學(xué)院,上海 201620; 2.上海市松江區(qū)新橋職業(yè)技術(shù)學(xué)校,上海 201612; 3.上海司南衛(wèi)星導(dǎo)航技術(shù)股份有限公司,上海 201801)

聲源定位系統(tǒng)[1-2]常用于設(shè)備的故障診斷[3],也可用于新產(chǎn)品的噪聲源定位,在語音人機(jī)交互中發(fā)揮著重要作用。聲源定位算法分為最大輸出功率波束形成法、高分辨率譜估計(jì)法和到達(dá)時間差(TDOA)的聲源定位方法[4]。前兩種方法多用于聲源方位角估計(jì),算法復(fù)雜度高且計(jì)算量大,而TDOA的聲源定位算法結(jié)構(gòu)簡單且計(jì)算速度快,常用于估計(jì)聲源三維坐標(biāo)。

TDOA分為時延估計(jì)和位置估計(jì)兩個步驟。第一步時延估計(jì)[5]的方法有廣義互相關(guān)方法[6]、自適應(yīng)最小均方方法[7]、自適應(yīng)特征值分解法[8]等。針對廣義互相關(guān)的時延估計(jì)算法,文獻(xiàn)[9]提出了一種改進(jìn)PHAT加權(quán)函數(shù)的方法,提升了時延估計(jì)性能。為了提升時延估計(jì)的抗噪能力,文獻(xiàn)[10]提出廣義互相關(guān)的SCOT/PHAT聯(lián)合加權(quán)時延估計(jì)算法。TDOA聲源定位算法的第二步是位置估計(jì)。在對聲源位置進(jìn)行解算的方法如下:在嵌入式的應(yīng)用中,為了追求計(jì)算簡單常采用直接計(jì)算的方法,但是其抗干擾能力和定位精度都不高。牛頓迭代法[11]也常被用于位置估計(jì),但其定位準(zhǔn)確性依賴于初始值,且運(yùn)算時間長。球形插值法[12]由于計(jì)算速度快,定位精度相對較高,得到了廣泛應(yīng)用。但是其抗噪性能較弱,不適用于對定位精度要求較高的場合。文獻(xiàn)[12]介紹了基于球形插值法求解聲源位置的TDOA方法,并將其應(yīng)用到了SRP-PHAT(Steered Response Power-Phase Transform)的優(yōu)化過程中,但并沒有針對TDOA算法進(jìn)行改進(jìn)。文獻(xiàn)[13]在五元三維立體麥克風(fēng)陣列的基礎(chǔ)上設(shè)計(jì)了一維空間搜索法計(jì)算聲源的位置,但是此法僅適用于三維立體的陣列。

本文為了解決平面麥克風(fēng)陣列的聲源三維坐標(biāo)定位問題,在廣義互相關(guān)算法作為時延估計(jì)算法的TDOA聲源定位算法中引入粒子群優(yōu)化算法,將其用于TDOA聲源定位的位置估計(jì)過程中,從而使TDOA聲源定位算法適用于平面麥克風(fēng)陣列。此外,本文研究發(fā)現(xiàn)相較于傳統(tǒng)的直接計(jì)算法、牛頓迭代法和球形插值法,基于粒子群優(yōu)化的TDOA聲源定位方法具有更好的魯棒性、定位精度和抗噪性能,能夠彌補(bǔ)TDOA聲源定位算法在位置解算過程中的不足。

1 時延估計(jì)

假設(shè)兩個麥克風(fēng)接收到的信號為xi(n)和xj(n)

xi(n)=s(n-τi)+ni(n)

(1)

xj(n)=s(n-τj)+nj(n)

(2)

其中,s(n)為聲源信號;ni(n)和nj(n)分別為兩個麥克風(fēng)接收到的噪聲;τi和τj分別為兩個麥克風(fēng)相對于參考麥克風(fēng)接收到的聲源時延,所需估計(jì)的時延差τij如式(3)所示。

τij=τj-τi

(3)

廣義互相關(guān)函數(shù)法[13]Rgcc(τ)為

(4)

其中,Gxixj=Xi(ω)Xj*(ω)為互功率譜密度函數(shù);Xi(ω)和Xj(ω)分別為xi(n)和xj(n)傅里葉變換后的結(jié)果;Xj*(ω)是Xj(ω)取共軛;ψij為加權(quán)函數(shù);e表示歐拉數(shù)。當(dāng)加權(quán)函數(shù)為PHAT函數(shù)時,如式(5)所示時,此加權(quán)方式在不改變相位信息的情況下,對功率譜函數(shù)進(jìn)行了白化處理,具有良好的抗噪性能。

(5)

求得Rgcc(τ)的值并對其進(jìn)行峰值檢測,即可得到時延估計(jì)結(jié)果。廣義互相關(guān)法求時延估計(jì)的原理如圖1所示。

圖1 廣義互相關(guān)法原理Figure 1.Principle of generalized cross-correlation method

2 粒子群優(yōu)化的位置估計(jì)算法

假設(shè)聲源在空間中的位置為Ss=[xs,ys,zs]T,第i個麥克風(fēng)的位置為mi=[xi,yi,zi]T,則聲源到達(dá)麥克風(fēng)j和麥克風(fēng)i的時間差τij表示為

(6)

式中,c表示聲速,一般取340 m·s-1。式(6)中聲源位置的坐標(biāo)是待求解的未知數(shù),至少需要3個方程聯(lián)列才能求得,以1號麥克風(fēng)為參考麥克風(fēng),選取4個麥克風(fēng),列出方程組。由于本文研究平面麥克風(fēng)陣列,所有麥克風(fēng)均分布在xoy平面,即zi=0,i=1,2,3,4。為了簡化計(jì)算,以(x1,y1,z1)作為聲源定位系統(tǒng)中三維笛卡爾坐標(biāo)系的原點(diǎn),列出如式(7)所示的方程組。

(7)

(8)

關(guān)于非線性多參數(shù)目標(biāo)函數(shù)優(yōu)化的問題,粒子群優(yōu)化方法相比于遺傳算法速度更快,相比于最速下降法更容易找到全局最優(yōu)解,因此本文引入粒子群優(yōu)化算法對聲源位置進(jìn)行搜索。

粒子群優(yōu)化[14-18]的步驟如下所示:

步驟1初始化粒子的位置和速度,多個粒子組成種群;

步驟2計(jì)算所有粒子的適應(yīng)度值F;

步驟3當(dāng)F小于粒子的歷史最佳位置對應(yīng)的適應(yīng)度時,用粒子的位置替換Pbest=(p1,p2,…,pd)中對應(yīng)粒子的最佳位置,Pbest為所有粒子的歷史最佳位置;

步驟4當(dāng)所有粒子的歷史最佳位置計(jì)算出的最小目標(biāo)函數(shù)值小于全局的歷史最佳目標(biāo)函數(shù)值時,用對應(yīng)最小目標(biāo)函數(shù)值的粒子位置替換全局歷史最優(yōu)位置Gbest=(g1,g2,…,gd);

步驟5按照式(9)更新粒子的速度,按照式(10)更新粒子的位置

(9)

xi(t+1)=xi(t)+vi(t+1)

(10)

其中,c1和c2分別為粒子學(xué)習(xí)因子和種群學(xué)習(xí)因子;r1和r2是兩個隨機(jī)數(shù);t是迭代次數(shù);vi=(vi1,vi2,…,vid)和xi=(xi1,xi2,…,xid)分別表示第i個粒子的速度和位置。在聲源定位的位置估計(jì)中,粒子的位置即聲源的位置,因此空間維數(shù)d為3。在式(9)中,慣性權(quán)重w的引入使粒子群算法搜索的過程發(fā)生了改變。w值越大,則vi(t)值越大,粒子在全局范圍內(nèi)搜索,適用于算法迭代初期,使得粒子能夠快速聚集到最優(yōu)解附近,加快算法的收斂;w值越小,則vi(t)值越小,此時粒子在局部范圍內(nèi)搜索,更擅長在小范圍內(nèi)接近最優(yōu)解,可提高算法收斂精度,適用于最后接近最優(yōu)解的迭代。因此,本文選用一種隨迭代次數(shù)線性減少的慣性權(quán)重的設(shè)置方式,其計(jì)算式為

(11)

式中,wmax和wmin分別為權(quán)重最大值和最小值;t和tmax為當(dāng)前迭代次數(shù)和最大迭代次數(shù);

步驟6滿足迭代停止條件,輸出全局歷史最優(yōu)位置,不滿足則返回步驟2。

基于粒子群優(yōu)化的聲源定位流程如圖2所示。

圖2 基于粒子群優(yōu)化的聲源定位的流程Figure 2.Flow of sound source localization based on particle swarm optimization

3 算法仿真分析

本文以如圖3所示的十六元平面陣列為聲源定位仿真的麥克風(fēng)陣列,以如圖4所示的汽車鳴笛的信號為聲源信號,對本文提出的基于粒子群空間搜索的聲源定位算法進(jìn)行仿真實(shí)驗(yàn)。

圖3 平面麥克風(fēng)陣列Figure 3. A planar microphone array

圖4 汽車鳴笛聲源信號Figure 4. Car whistle sound source signal

3.1 計(jì)算效率仿真

在處理器為Intel i7-9750H,內(nèi)存為32 GB的硬件條件下,分別對直接計(jì)算法、牛頓迭代法、粒子群空間搜索法以及球形插值法進(jìn)行了100次仿真實(shí)驗(yàn),各算法完成位置估計(jì)所用的平均時間如表1所示。

表1 各算法計(jì)算效率Table 1.Computational efficiency of each algorithm /s

在表1中,牛頓迭代法的迭代停止條件為誤差小于10-3時的計(jì)算時間。在基于粒子群優(yōu)化算法的空間搜索法中,將種群數(shù)量設(shè)置為1 000,迭代次數(shù)設(shè)為60,c1=c2=0.5,速度限制為[-1,1],慣性權(quán)重設(shè)為[0.40,0.95]。

由于直接計(jì)算法和球形插值法是根據(jù)理論推導(dǎo)出的計(jì)算式直接進(jìn)行求解,不涉及迭代過程,所以兩者速度最快。由表1可知,雖然引入了迭代過程,但是基于粒子群的空間搜索法與直接計(jì)算法、球形插值法計(jì)算速度相近。在牛頓迭代法中,初始值的選擇影響結(jié)果的準(zhǔn)確性,誤差精度的設(shè)置影響迭代次數(shù),對牛頓迭代法的計(jì)算速度產(chǎn)生影響,并且牛頓迭代法的計(jì)算涉及求導(dǎo)計(jì)算。相比于牛頓迭代法,基于粒子群的空間搜索法的迭代次數(shù)和種群數(shù)量可以根據(jù)實(shí)際要求進(jìn)行人為設(shè)置?；诹Ｗ尤旱目臻g搜索法計(jì)算效率相較于牛頓迭代法具有較大提高。

3.2 算法靈敏度分析

在TDOA的聲源定位算法中,聲源的定位結(jié)果受第一步時延估計(jì)的影響較大。在實(shí)際聲源定位中,時延估計(jì)會出現(xiàn)偏差,在安裝麥克風(fēng)時也會引入偏差。本文將兩者誤差統(tǒng)一轉(zhuǎn)換為各麥克風(fēng)與參考麥克風(fēng)聲程差的誤差。

(12)

(13)

δ表示所有M-1對麥克風(fēng)聲程差誤差的最大誤差,即在不同的最大誤差δ的情況下,所有麥克風(fēng)的誤差絕對值都要小于等于設(shè)定的最大誤差δ的絕對值。關(guān)于聲源定位的位置估算的結(jié)果誤差定義如下所示。

x軸坐標(biāo)的誤差為

(14)

y軸坐標(biāo)的誤差為

(15)

z軸坐標(biāo)的誤差為

(16)

聲源定位結(jié)果相對誤差為

(17)

隨著聲程差最大誤差δ的不同,聲源定位的結(jié)果誤差也隨之發(fā)生變化,本文部分實(shí)驗(yàn)設(shè)置最大誤差δ的范圍為[-2%,2%]。在此范圍內(nèi),麥克風(fēng)與聲源的距離偏差都在1 cm以內(nèi),以此來驗(yàn)證位置估計(jì)算法的靈敏度。

由于直接計(jì)算法較為簡單,性能較差,牛頓迭代法計(jì)算速度慢且依賴于初始值的選取,因此在以下仿真實(shí)驗(yàn)中不對這兩種方法進(jìn)行實(shí)驗(yàn)對比。球形插值法是基于最小二乘的閉式估計(jì)算法,可利用幾何關(guān)系計(jì)算的距離差和時延估計(jì)得到的距離差求得矩陣方程的閉式解。其誤差方程與式(8)形式類似,且球形插值法在TDOA聲源定位方法中應(yīng)用廣泛,因此在下述算法性能分析實(shí)驗(yàn)中,將粒子群空間搜索法與球形插值法進(jìn)行仿真實(shí)驗(yàn)對比。

圖5是在不同最大誤差下,分別利用粒子群空間搜索法與球形插值法對坐標(biāo)為(6 m, 5 m, 9 m)的聲源進(jìn)行位置估計(jì)所得到的位置誤差結(jié)果比較圖。

(a)

(b)

(c)

(d)圖5 位置估計(jì)誤差結(jié)果比較(a)x軸誤差 (b)y軸誤差 (c)z軸誤差 (d)相對誤差Figure 5. Comparison of position estimation error results(a)x-axis error (b)y-axis error (c)z-axis error (d)Relative error

由圖5可知,無論是粒子群空間搜索法還是球形插值法,兩者位置估計(jì)的誤差都隨最大誤差的增大而增大,其中粒子群空間搜索法在4種誤差中普遍小于球形插值法的誤差,并且大部分的誤差都能保證在3%的誤差范圍內(nèi),只有個別誤差會達(dá)到4%左右。而球形插值法隨著最大誤差的增大,其結(jié)果誤差變化較大,z軸的坐標(biāo)誤差甚至達(dá)到了10%左右。

以最大誤差為橫坐標(biāo),以位置估計(jì)的結(jié)果為縱坐標(biāo),得到如圖6所示的位置估計(jì)結(jié)果比較。由圖6(a)和圖6(b)可知,粒子群空間搜索的x軸和y軸定位結(jié)果和真實(shí)值的誤差在±0.2 m之間,而球形插值法則達(dá)到了-0.3 m的誤差。由圖6(c)可知,兩種算法z軸定位結(jié)果誤差比x軸和y軸坐標(biāo)的誤差更大,大部分誤差在±0.4 m的范圍內(nèi),但是球形插值的結(jié)果誤差會達(dá)到-0.6 m左右。綜合上述結(jié)果可知,粒子群空間搜索法的定位誤差在各方面都要小于球形插值法的定位誤差,具有較好的魯棒性。

(a)

(b)

(c)圖6 位置估計(jì)結(jié)果比較(a)x軸定位結(jié)果 (b)y軸定位結(jié)果 (c)z軸定位結(jié)果Figure 6. Comparison of position estimation results(a)x-axis positioning result (b)y-axis positioning result (c)z-axis positioning result

3.3 不同信噪比下的定位性能

將時延估計(jì)的計(jì)算結(jié)果加入到聲源定位的仿真實(shí)驗(yàn)中,在不同的信噪比下完成從時延估計(jì)到位置估計(jì)的整個聲源定位算法的仿真。隨機(jī)生成100個聲源位置坐標(biāo),其中x的取值范圍為[0 m,10 m],y的取值范圍為[0 m,10 m],z的取值范圍為[2 m, 20 m]。如圖4所示的汽車鳴笛信號為仿真信號,在不外加噪聲的情況下,粒子群空間搜索法和球形插值法的聲源定位誤差如圖7所示。

(a)

(b)

(c)

(d)圖7 仿真信號位置估計(jì)誤差結(jié)果比較(a)x軸誤差 (b)y軸誤差 (c)z軸誤差 (d)相對誤差Figure 7. Comparison of position estimation error results of simulation signal(a)x-axis error (b)y-axis error (c)z-axis error (d) Relative error

由圖7可知,粒子群空間搜索法的誤差在15%左右及以下,球形插值法的誤差大都超過了20%,表明粒子群空間搜索法的定位結(jié)果優(yōu)于球形插值法。將圖5中各誤差的均值和方差進(jìn)行統(tǒng)計(jì),結(jié)果如表2所示。

表2 無噪聲情況下的定位誤差Table 2. Positioning error without noise

由表2分析可得,粒子群空間搜索法的每一項(xiàng)誤差的均值都要比球形插值的誤差均值小5%左右。從方差部分可以看出,粒子群空間搜索法的誤差的穩(wěn)定性也優(yōu)于球形插值法。在無噪聲的情況下,兩種方法均出現(xiàn)了一定誤差,導(dǎo)致產(chǎn)生誤差的因素為時延估計(jì)算法在時延估計(jì)時的偏差和原始信號夾雜著噪聲。當(dāng)在仿真信號中人為地加入高斯噪聲的干擾,設(shè)置信噪比為5 dB、0 dB和-5 dB,粒子群空間搜索法和球形插值法的聲源定位誤差的均值和方差分別為表3～表5所示。

表3 信噪比為5 dB的定位誤差Table 3. Positioning error with signal-to-noise ratio of 5 dB

表4 信噪比為0 dB的定位誤差Table 4. Positioning error with signal-to-noise ratio of 0 dB

表5 信噪比為-5 dB的定位誤差Table 5. Positioning error with signal-to-noise ratio of -5 dB

由表3～表5可知,在信噪比從5 dB降到-5 dB的過程中,粒子群空間搜索法的各項(xiàng)誤差均值以大約10%的誤差進(jìn)行擴(kuò)大,而在此過程中,球形插值法的誤差均值隨著信噪比的降低而迅速增大。當(dāng)信噪比為5 dB時,球形插值法的誤差均值為40%～50%。當(dāng)信噪比為0 dB及以下時,球形插值法已經(jīng)無法有效地實(shí)現(xiàn)聲源定位。粒子群空間搜索法在-5 dB時的誤差均值要優(yōu)于球形插值法在5 dB時的誤差均值,表明粒子群空間搜索法在不同信噪比下的定位性能都要優(yōu)于球形插值法。

4 結(jié)束語

本文將粒子群優(yōu)化算法應(yīng)用于TDOA聲源定位算法的位置估計(jì)中,提出了以時延真實(shí)值和估計(jì)值差值的平方和為適應(yīng)度函數(shù)的粒子群算法的空間搜索法。

在計(jì)算效率方面,本文通過仿真對比了直接計(jì)算法、牛頓迭代法以及球形插值法,并在靈敏度和不同信噪比下的聲源定位性能這兩個方面對比了基于粒子群算法的空間搜索法和球形插值法。仿真結(jié)果表明,粒子群空間搜索法和球形插值法在計(jì)算速度上相近,且粒子群空間搜索法在魯棒性、準(zhǔn)確性和抗噪性能方面均優(yōu)于球形插值法。