基于遺傳算法的傳聲器陣列位置有源校正研究?

邵祚桪 厲 劍 鄭成詩 李曉東

(1 中國科學院聲學研究所 北京 100190)

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

0 引言

傳聲器陣列具有豐富的空時信息，能夠降低噪聲，抗干擾，進行靈活波束控制，實現(xiàn)聲源定位、識別、增強與分離等功能。傳聲器陣列按照幾何結構，可以分為一維線性陣列、二維平面陣列以及三維陣列，其中二維平面陣列廣泛應用于設備故障檢測、鳴笛檢測、遠距離拾聲等應用場景[1?2]。在應用過程中，傳聲器陣列的機械加工與安裝會導致位置不精確，產(chǎn)生平移與旋轉誤差等問題，使得傳聲器陣列相關算法在高精度定位應用場景達不到預期效果，因此需要研究陣列位置校正方法，從而提高定位性能。

針對傳聲器陣列位置校正問題，相關學者提出了一系列解決方法。文獻[3–5]使用已知位置的校正信號源，利用特定形狀傳聲器陣列的到達方向(Direction of arrival, DOA)信息，建立代價函數(shù)，通過迭代算法求解出各陣元的位置。文獻[6]利用陣列的到達時間(Time of arrival, TOA)建立模型，對陣列形狀和幅度進行校正。文獻[7–9]對到達時間模型的代價函數(shù)進行改進，并采用不同的改進優(yōu)化算法，提高了位置校正的準確性，并提高了收斂速度。上述基于TOA 方法需要將聲源與測量設備進行時鐘同步，實際應用中對于設備性能要求較高。文獻[10]通過建立最大似然方程，實現(xiàn)了陣元位置與聲源DOA信息的聯(lián)合估計。文獻[11]利用特征值分解法建立矩陣方程，利用最小二乘法，實現(xiàn)了陣元位置和幅相誤差的聯(lián)合估計，在高信噪比(Signal-noise ratio, SNR)下有較好的性能。文獻[12]利用多元變量分析對陣元幾何結構進行校正，并分析了算法對于陣列孔徑與位置的敏感度，驗證了較大的陣列尺寸可以保證校正準確性。在室內場景應用場景中，文獻[13]考慮晚期混響的影響，在多元變量分析的基礎上，引入多通道線性預測(Multi-channel linear prediction, MCLP)方法，提高了陣列位置估計的準確性。但是，上述方法在信號模型時均假設各陣元位置誤差彼此獨立，不能直接應用于形狀固定的陣列位置校正問題。

本文針對形狀固定的陣列位置校正問題，提出了一種基于遺傳算法的位置參數(shù)有源校正方法。首先，在校正聲源坐標已知的情況下，利用廣義互相關函數(shù)法(Generalized cross-correlation, GCC)估計每對傳聲器的到達時間差(Time difference of arrival, TDOA)并利用線性插值方法進行修正，構造均方誤差代價函數(shù)，再通過遺傳迭代算法最小化代價函數(shù)，得到全局最優(yōu)解。仿真及實驗結果表明，本方法校正誤差更小，且魯棒性較強，在陣列的各種應用場景具有普適性和實用性。

本文的主要內容如下：第1 節(jié)介紹了陣列與信號的模型；第2 節(jié)介紹了校正算法的具體原理與詳細步驟；第3 節(jié)給出陣列校正算法的仿真及結果分析；第4 節(jié)介紹了實驗及結果分析；第5 節(jié)為本文結論。

1 陣列模型介紹

考慮三維空間中一個由K個陣元組成的二維平面陣列，當存在J個時域或頻域可分的校正聲源時，第j個聲源向該陣列的第k個陣元入射，則該陣元接收到的信號頻域模型可表示為

其中，j=1,2,···,J，k=1,2,···,K，Sj(ω)表示第j個聲源的輸入信號，Hkj(ω)表示空間內第j個聲源到第k個陣元的傳遞函數(shù)，Wkj(ω)表示加性噪聲。

本文以陣列相對坐標對陣列傳聲器分布進行描述，以陣列所在平面為相對坐標yOz平面，以yOz平面內某點為坐標原點并由原點定義坐標軸，如圖1所示。設該陣列的相對坐標下第k個陣元的坐標為rin,k= [xin,k,yin,k,zin,k]，其中k= 1,2,···,K。以Tait-Byrne 角[14]的標準形式對該陣列進行旋轉變換，先后繞著相對坐標中的z軸、y軸、x軸旋轉角度α、β和γ，再將此陣列按照相對坐標中心在全局坐標中的位置平移rm= [x0,y0,z0]，則坐標變換后各陣元的全局坐標為

圖1 陣列的相對坐標模型Fig.1 The relative coordinate model of array

其中，

表示陣列分別繞z軸、y軸、x軸旋轉角度α、β和γ方式下的旋轉矩陣，其中，

因此描述該陣列在三維空間內平移與旋轉的待求未知參數(shù)集u可表示為

2 基于遺傳算法的陣列校正

使用迭代算法可對陣列參數(shù)集u各變量進行估計，本文提出采用遺傳算法[15]進行迭代求解，下面對校正算法的步驟進行闡述。

2.1 代價函數(shù)

遺傳算法需要設定代價函數(shù)表示適應度，針對本文的陣列校正問題，代價函數(shù)設為該陣列的每對陣元TDOA均方誤差函數(shù)。

本文陣列模型考慮陣列與校正源的三維坐標，因此將校正聲源視作球面波。由于陣列中各陣元位置不同，聲源到各陣元存在到達時間差。當?shù)趈個聲源信號入射時，對應的第l個陣元與第k個陣元(k≠l)接收信號分別為Xlj(ω)與Xkj(ω)，可利用GCC 函數(shù)估計兩路信號的到達時間差[16]。兩路信號的廣義互相關函數(shù)為

其中，ψ(ω)為GCC 的加權函數(shù)，Glkj(ω)為通道l和通道k信號的互功率譜。對廣義互相關函數(shù)進行峰值檢測可得到TDOA的估計值。由于使用GCC估計TDOA的精度受制于傳聲器的采樣率，為提高TDOA 估計的精度，可對互相關函數(shù)進行線性插值再進行峰值檢測。通過各陣元到聲源距離的幾何關系，可得到TDOA的理論值為

其中，rl和rk表示陣元l和陣元k的全局坐標，rj= [xj,yj,zj]表示第j個聲源的全局坐標，c為聲速。

將該陣列內每對陣元TDOA 的均方誤差相加，由此可得關于聲源j的均方誤差表達式：

再將所有J個聲源均方誤差相加，得到總體代價函數(shù)：

迭代過程中，代價函數(shù)越小，表明該個體適應度越高，該函數(shù)為遺傳算法篩選種群的唯一準則，代價函數(shù)對應的全局最優(yōu)解為

依據(jù)以上步驟，可得到迭代算法使用的代價函數(shù)，從而估計出各位置變量。

2.2 校正步驟

遺傳算法模擬生物學中的遺傳規(guī)律，對一個種群中的個體進行選擇、交叉、變異等過程進行概率化尋優(yōu)，求解全局最優(yōu)解?；谶z傳算法的陣列位置校正大致步驟如圖2所示，具體包括：

圖2 本文提出的陣列校正算法步驟Fig.2 The procedure of the proposed array calibration algorithm

(1)精確測量多個校正聲源的位置，使用陣列采集聲源信號，要求聲源時域或頻域可分。

(2)利用初步預計的陣列參數(shù)進行陣列坐標轉換，對陣列的TDOA 進行估計，作為遺傳算法的輸入變量。

(3)初步測量陣列各位置參數(shù)[x0,y0,z0,α, β,γ]，設置其中各變量的上下限，創(chuàng)建個體數(shù)為N的遺傳算法最初種群U=[u1,u2,···,uN]，種群中每個個體u1,u2,···,uN各位置參數(shù)由高斯分布產(chǎn)生，參數(shù)均值為初步測量所得數(shù)值，若超出預定上下限則重新生成，將每個個體的所有位置參數(shù)作為遺傳算法中的基因。設置遺傳算法的遺傳代數(shù)(及即迭代次數(shù))。

(4)對種群每個個體進行適應度(即代價函數(shù)Q)的計算。記錄種群最佳適應度、平均適應度與最佳的染色體。

(5)使用輪盤賭準則，對父代選出合適的染色體傳入下一代。

(6)使用單點交叉法，在基因中隨機尋找某個交叉點，使得某一對個體的基因進行彼此交換。

(7)對于個體的基因進行隨機的基本位變異。

(8)依據(jù)適應度選出上一次最佳與最劣適應度的染色體及其在種群的位置，新一次迭代最佳的染色體替代上一次最劣的位置，并記錄此時的平均適應度與最佳適應度。

(9)若達到最大遺傳代數(shù)，則停止迭代并輸出最佳染色體；若未達到則產(chǎn)生新的種群，編碼并進入步驟(4)繼續(xù)迭代，直至迭代結束。遺傳算法輸出結果即為陣列位置參數(shù)的全局最優(yōu)解。

3 仿真結果與分析

3.1 陣列位置校正方法比較

本文將上述方法與文獻[3]、文獻[17]中的方法進行比較，本文方法與其他文章方法均采用迭代算法求解目標變量，本文主要比較各方法代價函數(shù)的性能，優(yōu)化算法均使用2.2節(jié)所述的遺傳算法步驟。

文獻[3]為DOA-Based 方法，對特定形狀陣列接收到校正源j的信號，使用廣義互相關函數(shù)估計部分傳聲器對之間的TDOA，從而計算出相對坐標下的方位角與俯仰角，再利用陣列各陣元與聲源距離的幾何關系，計算出方位角與俯仰角θj與φj的理論值，其中j=1,2,···,J。每個聲源對應的角度均方誤差求和可得到代價函數(shù)：

文獻[17]為特征空間ES-Based(Eigenstructure-Based)方法，該方法利用多信號分類(Multiple signal classification, MUSIC)算法得到代價函數(shù)。設Rs,j(ω)為陣列各通道接收到校正源j的信號在頻率ω下的協(xié)方差矩陣，將協(xié)方差矩陣進行特征值分解，得到特征向量矩陣Uj= [uj,1,uj,2,···,uj,K]，其中uj,1為信號子空間的特征向量，uj,2,uj,3,···,uj,K為噪聲子空間的特征向量。由此得到噪聲子空間矩陣Un,j(ω)= [uj,2,uj,3,···, uj,K]。通過估計出的噪聲子空間矩陣?Un,j(ω)得到MUSIC 空間譜在聲源j對應的代價函數(shù)，將所有聲源對應的值相加得到總代價函數(shù)：

其中，Aj(ω)表示在頻率ω下陣列各陣元到聲源j的導向矢量，符號“(·)H”表示矩陣的共軛轉置，符號“‖·‖”表示矩陣的F-范數(shù)。

為比較本文方法與上述文獻方法的性能，本文選取二維平面?zhèn)髀暺麝嚵械膽脠鼍爸坏镍Q笛聲定位場景[18]進行校正仿真，使用4×8 的平面均勻矩形陣列，陣元共32個，相鄰陣元相距0.04 m，將該陣列角上的陣元位置定義為該陣列的相對坐標原點。

校正源使用信號為時長2 s 的鳴笛信號，對應語譜圖如圖3所示。本文仿真中在各傳聲器加入互不相關的白噪聲，測試不同信噪比下的各種算法的校正性能。本文考慮全頻帶定位方法下的陣列校正問題，在仿真測試中公式(8)中ψ(ω)在所有頻率上均取1。

圖3 鳴笛的語譜圖Fig.3 The spectrogram of the whistle

鳴笛的戶外場景應用以環(huán)境噪聲為主而混響較低，本文在仿真中只考慮直達聲模型，即公式(1)中的傳遞函數(shù)Hkj(ω)取1。面陣原始位置(即上述陣元的相對坐標原點的位置)[x0,y0,z0,α,β,γ]參數(shù)分別為(1 m, 1 m, 8 m, 0°, 0°, 0°)，坐標參數(shù)搜索范圍為原位置±0.1 m，角度參數(shù)搜索范圍為原位置±5°。聲源坐標分別為(10 m, 5 m, 1 m)、(10 m,0 m, 1 m)、(20 m, 5 m, 1 m)，陣列與校正源位置的示意圖如圖4所示。按照全局信噪比30 dB、20 dB、10 dB 三種情況仿真，每種信噪比情況進行20 次Monte-Carlo測試，其余仿真參數(shù)如表1所示。

表1 仿真參數(shù)Table 1 Simulation parameters

圖4 仿真場景示意圖(單位：m)Fig.4 The schematic diagram of the simulation scene(unit: m)

圖5給出了本文方法(Proposed)、DOA-Based方法、ES-Based 方法對于陣列坐標參數(shù)與陣列旋轉角度參數(shù)的平均校正誤差。圖5(a)～(c)結果表明，在信噪比10 dB 及以上的多數(shù)情況下，對坐標參數(shù)的估計，本文方法誤差小于DOA-Based算法與ES-Based 算法，最大誤差限制在0.016 m 以內。圖5(d)～(f)結果表明，在各信噪比下，本文方法對角度參數(shù)的估計誤差也總體小于其他兩種算法，最大誤差限制在0.08°以內。

圖5 不同信噪比下各方法陣列參數(shù)校正平均誤差Fig.5 The average calibration errors of array parameters with different SNRs

經(jīng)分析，DOA-Based 算法使用遠場模型，僅考慮角度變量，未考慮各陣元至聲源的距離變量，影響了代價函數(shù)的準確性。ES-Based 方法中MUSIC 算法空間譜在特定方向上可形成極高的譜峰，但該方法對噪聲較敏感，噪聲環(huán)境下聲源方向空間譜的實際值與理論值存在偏差，影響迭代中全局最優(yōu)解的搜索。本文方法采用近場模型，考慮陣列與聲源的三維坐標，同時在GCC 函數(shù)中引入線性插值，提高了低信噪比下TDOA估計的準確性，使得迭代過程中陣列各位置參數(shù)更容易收斂至預期位置。

綜上所述，相比于傳統(tǒng)的DOA-Based 和ESBased 方法，本文提出的方法可以獲得更加準確的陣列位置參數(shù)的估計性能。

3.2 陣列存在不一致時的校正性能

傳聲器陣列在制造或使用的過程中會產(chǎn)傳聲器之間幅度響應或相位響應不一致的情況，為驗證傳聲器不一致情況下本文算法的魯棒性，本文進行如下對比仿真，其中陣元不一致的情況包括：

(1)假設各個傳聲器陣元的幅度存在?1～1 dB均勻分布的誤差；

(2)假設各個傳聲器陣元的相位存在?3° ～3°均勻分布的誤差；

(3)假設各個傳聲器陣元的相位存在?3° ～3°均勻分布的誤差，且幅度存在?1～1 dB 均勻分布的誤差。

仿真中其他參數(shù)條件、遺傳算法參數(shù)均與3.1節(jié)所述相同。表2～4 分別給出了情況(1)、(2)、(3)對應的位置參數(shù)校正平均誤差，其中誤差最小的數(shù)據(jù)進行了加粗處理便于顯示。由仿真結果可知，與其他兩種算法相比，本文提出的校正算法在絕大多數(shù)情況下的參數(shù)誤差均為最小，因此本文算法在傳聲器陣列幅度與相位不一致的情況下具有更好的魯棒性。

表2 傳聲器陣列幅度不一致時參數(shù)校正平均誤差Table 2 The average calibration errors of array parameters when the amplitudes of different microphones are inconsistent

3.3 聲源位置對校正性能的影響

在實際應用中，聲源位置的選取會影響到傳聲器陣列位置校正的準確性，為研究聲源位置本文所提出校正方法性能的影響，本節(jié)采取不同的聲源位置組合分別進行仿真。

本節(jié)仿真假設仿真場景模型為直達聲模型。陣列原始位置(1 m, 1 m, 8 m, 0°, 0°, 0°)，坐標參數(shù)搜索范圍為原位置±0.1 m，角度參數(shù)搜索范圍為原位置±5°，仿真參數(shù)與表1相同。仿真在50 m×20 m的矩形區(qū)域進行3 個聲源的放置，這里將x軸方向定義為主要的距離方向，x軸數(shù)值越大代表距離陣列越遠，位置組合包括：

(A)近距離集中布放：(1 m, 1 m, 1 m)，(1 m,2 m, 1 m)，(2 m, 1 m, 1 m)；

(B)遠距離分散布放：(50 m,10 m,1 m)，(50 m,15 m, 1 m)，(50 m, 20 m, 1 m)；

(C)遠距離集中布放：(50 m, 1 m, 1 m)，(50 m,2 m, 1 m)，(50 m, 3 m, 1 m)；

(D)中距離分散布放：(20 m,10 m,1 m)，(25 m,10 m, 1 m)，(30 m, 10 m, 1 m)；

(E)不同距離分散布放：(1 m,1 m,1 m)，(20 m,10 m, 1 m)，(50 m, 20 m, 1 m)。

表5給出了不同聲源位置組合下各參數(shù)的平均校正誤差，在6 個位置參數(shù)在3 種信噪比條件下共有18 種情況。由仿真結果可知，校正聲源位置的選擇會影響各個參數(shù)的估計誤差結果。對比5 種布放組合的18個校正結果可知，組合(E)在11個變量上平均誤差小于其他組合，且多數(shù)情況為旋轉角度校正誤差較小。這是因為聲源位置分布較為分散的情況下，陣列微弱的角度偏轉對TDOA 的影響較大，使得以TDOA 均方誤差作為代價函數(shù)的校正算法對旋轉角度的估計更準確。此外，當聲源布放在x軸方向較遠情況時(即組合(B)和組合(C))，平移坐標的估計誤差總體大于其他情況，這是因為聲源距離陣列較遠的情況下，陣列TDOA值的變化量對陣列的平移不敏感，導致錯誤陣列位置對應的代價函數(shù)與正確位置相差較小，因此迭代算法容易收斂錯誤。

表5 不同聲源位置組合下的參數(shù)校正平均誤差Table 5 The average calibration errors of array parameters when using different combinations of sound sources with different positions

基于上述仿真結果，在聲源個數(shù)有限的情況下使用本文方法進行陣列校正時，應采用聲源位置較為分散布放方式提高參數(shù)校正精度。

3.4 其他形狀陣列校正

為驗證本文方法適用于不同形狀陣列的位置校正，本文另選取3 種形狀的陣列進行校正算法仿真，陣列形狀包括：

(a)十字陣列，陣元共32個，每一分支上相鄰陣元間距0.04 m；

(b)螺旋陣列，陣元共32個，螺旋每層半徑相距0.04 m；

(c)隨機陣列，陣元在坐標范圍x ∈[?0.2,0.2]、y ∈[?0.2,0.2]內均勻分布產(chǎn)生，陣元共32個。

圖6給出了3 種陣列在相對坐標下各陣元的位置。

圖6 不同陣型示意圖Fig.6 The schematic diagrams of arrays with different formations

仿真條件與3.1 節(jié)相同，假設仿真場景模型為直達聲模型。陣列原始位置(1 m, 1 m, 8 m, 0°, 0°,0°)，坐標參數(shù)搜索范圍為原位置±0.1 m，角度參數(shù)搜索范圍為原位置±5°。聲源位置(10 m, 5 m,1 m)，(10 m, 0 m, 1 m)，(20 m, 5 m, 1 m)。信噪比按照30 dB、20 dB、10 dB 三種情況仿真，每種信噪比情況進行20 次Monte-Carlo 測試，其余仿真參數(shù)如表6所示。

表6 仿真參數(shù)Table 6 Simulation parameters

表7給出了十字陣列(a)、螺旋陣列(b)、隨機陣列(c)的平均校正誤差。仿真結果表明，在10 dB及以上的信噪比情況下，本文方法在多數(shù)物理量估計的誤差都小于DOA-Based方法與ES-Based方法，本文提出的方法適用于多種類似陣列，且陣列位置的估計性能優(yōu)于其他兩種方法。

同時由表7可知，在使用本文方法情況下，坐標參數(shù)平均誤差滿足a

表7 不同信噪比下各形狀陣列參數(shù)校正平均誤差Table 7 The average calibration errors of array parameters for the arrays of different formations with different SNRs

4 實驗驗證

本節(jié)利用矩形陣列分析上述各個陣列位置校正算法在半消聲室中的性能。測試所使用的半消聲室長寬高尺寸為14.8 m×8.8 m×7.2 m。實驗使用的陣列為4×8的矩形均勻面陣，陣元間距為0.04 m，采樣率為16 kHz。錄聲使用專用系統(tǒng)綜合測試軟件控制采集板。校正聲源為惠威H4 揚聲器播放的鳴笛聲，信號源與圖3所示相同。陣列中心位置(2 m,4 m, 1.33 m)，校正源位置(2 m, 4 m, 1.12 m)，(3 m,4 m, 1.12 m)，(4 m, 4 m, 1.12 m)，聲源與采集板擺放位置的示意圖如圖7所示，場景實際拍攝圖片如圖8所示。

圖7 實驗場景示意圖(單位：m)Fig.7 The schematic diagram of the experiment scene(unit: m)

圖8 實驗場景圖片F(xiàn)ig.8 The photo of the experiment scene

實驗步驟如下：

(1)建立空間直角坐標系，確定坐標軸方向，設置3個聲源，準確測量聲源位置；

(2)擺放陣列，使得采集板大致與y軸方向垂直，大致測量其位置，定義此時其繞著z軸的偏轉角度為0°；

(3)在預定聲源位置先后分別播放2 s 鳴笛，使用采集板進行錄聲；

(4)將陣列繞z軸分別旋轉5°、10°、15°，重復步驟(3)。

利用4種陣列角度下陣列分別接收到的校正源信號，估計出每次的α角度，以第一輪實驗位置為基準0°，后3 次陣列旋轉情況下估計出的α與原始位置的α做差，得到Δα，每次角度重復估計5次，求出角度估計的均值與標準差。實驗所使用參數(shù)如表8所示。

表8 實驗參數(shù)Table 8 Experiment parameters

該實驗對本文提出方法、文獻[3]所述DOABased 方法及文獻[17]所述ES-Based 方法的性能進行測試。圖9的誤差棒圖給出了本次實驗中3 種校正方法對Δα的估計結果，其中點線與符號“o”代表真實的旋轉角度(True)，實線代表本文方法下旋轉角度估計的平均值與標準差，虛線代表DOABased方法估計結果，點劃線代表ES-Based 方法估計結果。

圖9 陣列旋轉角度估計結果Fig.9 The estimated results of array rotations for different angles

由實驗結果可知，在不同旋轉情況下的角度估計中，本文方法估計誤差總體低于其他方法，最大誤差僅0.76°，具有較高的準確性。此外，本實驗在不同α角度上多次重復估計下的標準差也明顯小于其他兩種方法，最大標準差僅為0.32°，這說明利用遺傳算法進行校正，在相同參數(shù)與迭代次數(shù)的情況下，本文方法所使用代價函數(shù)更易收斂至全局最優(yōu)解。綜上，該實驗的結果證明了本文方法具有更好的校正性能。

5 結論

本文提出了一種對二維平面?zhèn)髀暺麝嚵械淖鴺藚?shù)與旋轉角度參數(shù)校正的方法。首先利用TDOA估計值與理論值構造均方誤差代價函數(shù)，再使用遺傳算法搜索全局最優(yōu)解使得代價函數(shù)最小，由此求解陣列的位置參數(shù)。本文在不同信噪比情況下，對不同形狀的陣列進行了位置參數(shù)估計的仿真，并進行半消聲室的實測數(shù)據(jù)驗證。仿真與實驗結果表明，使用該代價函數(shù)并使用遺傳算法進行優(yōu)化可得到較為準確的位置參數(shù)，該方法在不同信噪比情況下校正誤差總體小于DOA-Based 算法與ES-Based 算法，并且適用于不同形狀的傳聲器陣列，收斂結果穩(wěn)定，魯棒性更強。