高 玥 盧 鈴 吳 鳴 楊 軍 曹 寅

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

(2 中國科學院噪聲與振動重點實驗室(聲學研究所)北京 100190)

(3 國網(wǎng)湖南省電力有限公司電力科學研究院 長沙 410007)

(4 國網(wǎng)電力設施噪聲與振動實驗室 長沙 410007)

(5 機器視覺 語音及信號處理研究中心 吉爾福德 GU2 7XH)

0 引言

傳聲器陣型的配置在解決聲成像問題中起著重要的作用。傳聲器陣列的設計在很大程度上取決于所期望的應用。Meyer 等[1]在2002年提出基于球諧分解的剛性球陣波束形成器，可在不改變波束模式的情況下將觀測方向指向三維任意位置。Rafaely[2]提出了一種基于球諧的球陣列設計和分析框架，還提出了球面上傳聲器定位的備選空間采樣方案。Balmages等[3]提出了開放雙球面陣，解決了由于貝塞爾函數(shù)零點引起的空間模態(tài)對應頻率上出現(xiàn)的數(shù)值病態(tài)問題。Carious 等[4]和Lamotte等[5]利用了雙球面陣列，其內部為剛性陣列，非常適合中高頻范圍，而外部是一個開放球面陣列，目的是提高低頻性能，平衡了內部陣列方向性好但分辨率差和外部陣列的分辨率好但方向性差的問題，從而實現(xiàn)較好性能。上述研究在低頻段均存在著所需陣列孔徑較大不易制作和運輸?shù)膯栴}。

基于傳聲器陣列的聲成像方法以其簡單和高效的優(yōu)勢有著廣闊的應用背景，可分為陣元域和模態(tài)域兩種。傳統(tǒng)陣元域的算法中一類為間接法，其利用不同陣元接收信號的到達時間差來估計聲源方位。其中，廣義互相關相位變換方法使用廣義互相關方法進行時間延遲估計[6]，并對幅度譜引入加權函數(shù)做相位變換[7]進行聲源方位估計，由于其復雜度較低因此得到了廣泛的發(fā)展和應用，但其受混響的影響較大。在此基礎上，Michaud 等[8]研究了實時追蹤移動聲源的問題，Grondin 等[9]在噪聲條件下通過時頻掩膜技術來提高定位精度。另一類直接法是計算所有方位上的代價函數(shù)從而估計聲源位置。其中可控響應功率法[10]的具體實現(xiàn)與波束形成器的設計方法相關，在智能音箱等領域有著較大應用，研究重點在于減少空間掃描的網(wǎng)絡格點數(shù)；最大似然法[11]處理相干信號能力較好且統(tǒng)計特性優(yōu)秀但是計算復雜度高；基于子空間的方法擁有超分辨率，其中非相干子空間法[12]的空間相干矩陣缺秩問題會導致算法性能下降，而相干子空間法則[13]需要一定的先驗知識來獲得聚焦矩陣。

傳統(tǒng)的波束形成算法在低頻段效果較差，且在相干聲源、復雜聲源存在的條件下算法性能會受限，因此本文采用廣義逆波束形成(Generalized inverse beamforming, GIBF)算法來解決上述問題[14?15]。廣義逆波束形成的改進主要以改進正則化方法為特征，求解反問題所涉及的優(yōu)化參數(shù)對于準確的波束形成結果至關重要，Zavala 等[16]對此進行了討論。為了驗證算法并評估其性能，Bahr 等[17]將廣義逆波束形成應用于與小型露天噴氣式設施即NASA蘭利靜音流動設備相對應的實驗基準數(shù)據(jù)集中，研究了由不同貢獻者實施并應用于同一測試用例的幾種陣列分析方法的可變性。與其他常見波束成形技術的比較可以顯示廣義逆波束形成在復雜應用中的潛力和優(yōu)勢。

在上述陣元域方法的基礎上，模態(tài)域方法根據(jù)傳聲器陣列的拓撲結構，把接收到的信號從陣元域進一步展開到模態(tài)域，將環(huán)形陣列信號轉換到環(huán)諧域，球形陣列信號轉換到球諧域。由于模態(tài)域波束形成器的導向矢量具有頻率無關的特性，可以直接實現(xiàn)頻域平滑從而減少計算量，對于實際應用有很大的價值。Rafaely[18]在球諧域采用延時求和的方法設計波束形成器，Gao 等[19]將球諧域聲源定位方法應用在水下環(huán)境中并分析其性能。

本文的核心思路為：提出了基于球諧分解的L2范數(shù)廣義逆波束形成算法，并在分布式開放球陣的陣型基礎上進行仿真及實驗。得到如下結論：通過仿真及實驗分析對比了不同算法在低頻相干聲源條件下聲成像的準確性和對陣元位置誤差的魯棒性后，證明在此分布式開放球陣陣型下該算法具有較好的性能。

1 聲成像算法原理

1.1 球諧系數(shù)估計法

Ueno 等[20]提出了適用于任意陣型和陣元指向性的全局諧波系數(shù)估計方法。通過貝葉斯準則，將截斷后的后驗均值作為區(qū)域中球諧波系數(shù)矢量的估計：

其中，區(qū)域中的任意一點用r=(r,θ,φ)表示，(r0,k)為在r0處展開的截斷后的球諧波系數(shù)矩陣，k為波數(shù)，P為截斷階數(shù)；λ是與球諧波系數(shù)矢量的先驗分布有關的常數(shù)，xe為接收的三維聲場，Σ ∈CQ×Q是噪聲的互功率譜密度矩陣。陣列由Q個陣元構成，分別分布在r1,···,rQ處，陣元指向性為c1,···,cQ，c是由球諧波分解系數(shù)cmn構成的矢量。ΞP(r0,k)是與不同展開中心展開球諧波系數(shù)之間的轉換矩陣T及陣元指向性系數(shù)c(θ,φ)有關的矩陣，Ψ(k)可由ΞP(r0,k)表示[20]，即：

其中，(·)?表示復共軛轉置。

考慮幅度為S(w)的單頻平面波從(θ0,φ0)方向入射，則在r處的頻域觀測信號的球諧波展開為[21]

其中，波矢k=(k,θ0,φ0)，jn(kr)為n階球貝塞爾函數(shù)，Ymn(θ,φ)為球諧波函數(shù)。由于任意三維聲場x(r,k)可以表示為球諧波展開的形式：

其中，αmn為球諧波系數(shù)，基函數(shù)定義為φmn(r,k)=則由公式(4)、公式(5)可得到若是D個平面波入射時，聲場中的球諧波系數(shù)具有如下形式：

其中，Sd(ω)和(θd,φd)分別為dth 平面波的幅度和入射方向。

根據(jù)球諧波函數(shù)的正交性， 定義濾波器h(θ,φ)∈C∞，令濾波器系數(shù)為

則該濾波器輸出為

可以看出，該濾波器具有理想的空間選擇特性。因此通過改變?yōu)V波器的掃描角度則可得到方位譜的估計值，從而確定信號的入射角度。

1.2 廣義逆波束形成

廣義逆波束形成可以歸類為逆方法，因為它需要對聲源-傳聲器間的傳播問題進行反演，可以建模為[22]

其中，A ∈CQ×N，q ∈CN×1，p ∈CQ×1，N為掃描網(wǎng)絡格點數(shù)。由于一般情況下Q ?N，則問題轉化為求欠定方程的解。本文采用一范數(shù)及二范數(shù)兩種方法來求解。

1.2.1 L1范數(shù)廣義逆波束形成

L1范數(shù)的廣義逆波束形成(L1-Generalized inverse beamforming, L1-GIBF)算法是用壓縮感知(Compressive sensing, CS)來對稀疏的噪聲源進行聲成像的算法[23]。CS 為求解欠定問題Aq=p提供了一種解決方案，其中測量矩陣p、聲源振幅q和線性傳遞矩陣A(即CS 中的感知矩陣)分別為Q×1 的向量、N ×1 的向量以及Q×N維矩陣，其中Q為組成傳聲器陣列的陣元數(shù)，N為格點數(shù)，且一般情況下Q ?N。陣元處聲信號可以由向量和矩陣形式表示為

其中，ε是由電噪聲、實驗誤差、測量數(shù)據(jù)與模擬聲壓不匹配等引起的誤差。A的列向量an等價于球面波束形成的導向矢量。由于Q遠小于N，導致求解x的欠定問題沒有唯一解。由于聲源位于特定的局部區(qū)域內，因此解決不適定問題的一種可行方法是將稀疏度放在x上。則可將求解式(10)轉化為約束問題：

但是L0范數(shù)的最小化問題是一個非凸問題，在計算上難度較大。當有足夠稀疏的聲源時，公式(11)等價于L1范數(shù)的最小化問題，則可以轉化為凸優(yōu)化問題從而有效地實現(xiàn)計算處理。又由于q中的一個元素乘以具有較大歐氏范數(shù)的A中的一列，在計算L1 范數(shù)最小化問題時趨向于一個非零解。為了防止這種偏差，將陣元處測量數(shù)據(jù)p和an分別按照歐幾里得范數(shù)進行歸一化?；谝陨峡傻玫?/p>

借助數(shù)學工具中凸優(yōu)化算法建模軟件包對式(12)進行求解。由于它的稀疏性，即使在單頻情況下，此算法也能提供高分辨率的聲成像結果。然而，當測量數(shù)據(jù)受到噪聲污染時，稀疏解的非零分量會出現(xiàn)在與實際源位置不同的位置。

1.2.2 L2范數(shù)廣義逆波束形成

L2 范數(shù)的廣義逆波束形成算法(L2-Generalized inverse beamforming, L2-GIBF)是用Tikhonov 正則化[24]修正問題的方法來求解公式(9)欠定方程的解。通過添加殘差范數(shù)項和正則化參數(shù)λ2修正等式后，得到正則化的解為

其中，A=USV?為奇異值分解結果，U ∈CQ×Q是酉矩陣，S ∈RQ×Q是一個對角矩陣，其中包含的奇異值按降序s1≥s2≥···≥sQ≥0 排列，V ∈CN×Q滿足V?V=IQ，正則化系數(shù)λ2定義了濾波器的強度。

1.3 基于球諧分解的L2范數(shù)廣義逆波束形成

基于球諧分解的L2 范數(shù)廣義逆波束形成(Spherical harmonics decomposition-L2-Generalized inverse beamforming, SHD-L2-GIBF)算法的原理是先通過球諧分解(SHD)的預處理來穩(wěn)定后續(xù)的求逆數(shù)值問題，尤其在中低頻范圍內效果更好，而后通過廣義逆波束形成(GIBF)處理截斷后的陣元接收數(shù)據(jù)進行聲成像。

1.3.1 算法原理

pnm為復數(shù)球諧波系數(shù)，是由球傅里葉變換(Spherical Fourier transform, SFT)計算得到的球諧函數(shù)的譜，表示為

其中，p(θq,φq)為第q個陣元的聲壓分布，αnmq是由積分轉換到累加和形式的近似，對第q個陣元的m階n次項的加權系數(shù)。由于球形陣列傳聲器個數(shù)Q是有限的，因此p(θq,φq)可以通過截斷P階后的逆球傅里葉變換(Inverse spherical Fourier transform,ISFT)近似表示為[25]

則可將式(15)轉換成矩陣形式表示為

其中，p為各個陣元采集聲壓數(shù)據(jù)經(jīng)過短時傅里葉變換的頻譜，pnm為球諧域譜，YP矩陣的階數(shù)為M ×(P+1)2，包含在每個陣元上計算得到P階的球諧函數(shù)Ymn(θq,φq)。

對式(16)中的線性系統(tǒng)進行反演，得到了SFT的數(shù)值解為式(17)：

其中，TP為YP的逆矩陣。

短時傅里葉變換后得到的信號頻譜先經(jīng)過ISFT 轉換到模態(tài)域，然后通過SFT 轉換回來，則實現(xiàn)了球諧分解步驟[25]。但是SFT和ISFT只能用于單球陣陣列條件下，為了克服這一限制，本文提出將球諧系數(shù)估計法引入此步球諧分解操作中。因此公式(16)、公式(17)可以用公式(1)、公式(8)的形式表示，則可在任意陣型下實現(xiàn)球諧分解預處理步驟。將預處理后的數(shù)據(jù)再進行本文1.2.2 節(jié)中敘述的二范數(shù)廣義逆波束形成算法，實現(xiàn)基于球諧分解的L2范數(shù)廣義逆波束形成算法。

1.3.2 正則化參數(shù)選擇

對于公式(13)中正則化參數(shù)的選取方法，本文采用的為L-曲線準則[26]，在二維上，以lg-lg 為尺度，(‖p ?Aq‖,‖q‖)的值構成了如圖1所示形狀類似于L 的曲線，故被稱為L 曲線。為了平衡欠正則化和過正則化，在L曲線曲率最大的點(即由豎直到水平的轉角)處選擇正則化參數(shù)。以λ2為參數(shù)的L曲線的曲率函數(shù)可以表示為

圖1 L 曲線圖Fig.1 L–curve

其中，u(λ2)= lg‖p ?Aq‖，v(λ2)= lg‖q‖，且“′”，“′′”表示對λ2求解一階導數(shù)和二階導數(shù)。

2 聲成像算法性能仿真分析

2.1 仿真陣型及條件

傳聲器陣列陣型的配置在解決聲成像問題中起著重要的作用。本文的研究重點在于低頻相干聲源的聲成像性能。對于固定尺寸的球形陣列，其定位極限也隨之確定，這種限制尤其在低頻段影響更為明顯。增加陣列尺寸能夠提高陣列的空間分辨率，但是在實際操作中隨著頻率降低不斷增大陣列孔徑是不現(xiàn)實的，而且大孔徑球形陣列的系統(tǒng)過于龐大，在制作、運輸和使用等方面都有很大困難。本文采用了一種分布式的開放球形陣列，子陣孔徑較小適用于中高頻范圍，外部整體陣列孔徑較大可以擴大頻率范圍到低頻。具體的布放方案為：為實現(xiàn)多聲源聲成像的目的，選取截斷階數(shù)P為4。傳聲器的分布方案要求符合奈奎斯特采樣定理的空間球諧形式Q≥(P+1)2，從而保證空間采樣不混疊，因此選取傳聲器數(shù)量Q= 32。32 陣元整體分布式球形陣列示意圖見圖2(a)，它由8 個子陣組成，子陣的球心位于正六面體的8 個頂點上，正六面體的外接球半徑為R= 0.6 m。子陣如圖2(b)所示，陣元位于正四面體的4 個頂點上，正四面體的外接球半徑為r=0.13 m。

圖2 陣列結構示意圖Fig.2 Geometry model of the array

仿真陣型為上述分布式開放球形陣列，仿真條件為空間中兩個等能量白噪聲聲源分別位于(θ1=?135°,φ1= 75°)和(θ2= 60°,φ2= 120°)，其中θ為方位角，φ為俯仰角；空氣中聲速c= 343 m/s，采樣率為16 kHz，幀長為640，在時域上添加一定信噪比(Signal-noise ratio, SNR)的高斯白噪聲為仿真噪聲。格點采集于方位角?180° ～180°、俯仰角0° ～180°的區(qū)間內，步長均為3°。仿真對比研究了6 種聲成像算法，包括常規(guī)頻域波束形成(Conventional frequency domain beamforming, CFDBF)[27]、函數(shù)波束形成(Functional beamforming, FUNBF)[28]、增強超分辨率空間聲源相干CLEAN(Enhanced high resolution-CLEAN-source coherence, EHR-CLEAN-SC)[29]算法、L1-GIBF、球諧系數(shù)估計法以及SHD-L2-GIBF。

2.2 相干聲源聲成像

兩個等能量相干白噪聲聲源， 在信噪比15 dB 條件下，300 Hz 處的6 種算法譜圖結果見圖3?？梢娫趦上喔陕曉创嬖诘臈l件下，低頻段處理時CFDBF 及以CFDBF 為基礎的兩種算法即FUNBF 和EHR-CLEAN-SC 的譜圖具有一定的相似性。CFDBF 由于其波束主瓣寬導致聲成像結果有虛峰且譜峰位置不準確。而FUNBF 雖然經(jīng)過參數(shù)冪次處理后旁瓣降低，但仍然存留著CFDBF 的譜峰不準確問題。EHR-CLEAN-SC 算法是基于旁瓣與主瓣相干的假設從而降低旁瓣，因此在低頻段相干聲源存在時性能大大降低。

圖3 在SNR=15 dB、300 Hz 仿真條件下相干聲源聲成像歸一化譜圖Fig.3 Normalized spectrogram of coherent sound source acoustic imaging under SNR = 15 dB,300 Hz simulation conditions

在SNR = 15 dB、300 Hz 條件下，通過50 次Monte-Carlo實驗，6 種算法對兩相干聲源定位誤差的均值結果見表1。

表1 SNR=15 dB、300 Hz 下，6 種算法對兩相干聲源50 次仿真定位誤差的均值Table 1 The mean of the localization error of the six algorithms for 50 times simulation of two coherent sound sources at f =300 Hz, SNR = 15 dB

對于低頻相干可成像的3 種算法進行進一步計算誤差方差可得到表2。

表2 SNR=15 dB、300 Hz 下，3 種算法對兩相干聲源50 次仿真定位誤差的方差Table 2 The variance of the localization error of the three algorithms for 50 times simulation of two coherent sound sources at f =300 Hz, SNR=15 dB

2.3 陣元位置誤差

由于在實際應用中，陣元位置誤差可能會引起一定程度上的算法性能損失。因此仿真引入每個陣元存在均值0、方差為1 cm的隨機位置誤差，再進行上述3 種低頻相干聲源聲成像算法50 次仿真計算可以得到定位誤差均值和方差，如表3和表4所示。

表3 SNR=15 dB、300 Hz 下，在隨機陣元位置誤差下，對兩相干聲源50 次仿真定位誤差的均值Table 3 The mean of the localization error of the three algorithms for 50 times simulation of two coherent sound sources at f =300 Hz,SNR=15 dB with random microphone position error

表4 SNR = 15 dB、300 Hz 下，在隨機陣元位置誤差下，對兩相干聲源50 次仿真定位誤差的方差Table 4 The variance of the localization error of the three algorithms for 50 times simulation of two coherent sound sources at f = 300 Hz, SNR = 15 dB with random microphone position error

通過上述分析得到結論，在此分布式開放球形陣列陣型配置下，低頻相干聲源時CFDBF、FUNBF、EHR-CLEAN-SC算法由于在低頻導向矢量的相干性較強，波束的主瓣較寬，導致聲源位置估計誤差較大，無法達到聲成像要求。SHD-L2-GIBF、L1-GIBF 和球諧系數(shù)估計法有較好的性能。其中L1-GIBF 算法由于要進行最優(yōu)化迭代，擁有較高分辨率，可用于區(qū)分距離較近的聲源，但是計算復雜度較高，且算法穩(wěn)定性較差。球諧系數(shù)估計法由于模態(tài)展開的頻率無關性，可以直接實現(xiàn)頻域聚焦，擁有較好的定位性能。SHD-L2-GIBF 的聲成像由于球諧分解預處理的截斷作用，高頻的誤差及噪聲對算法的影響效果減弱，為算法提供了數(shù)值穩(wěn)定性。通過多次實驗得到誤差均值及方差的表格對比，可以定量得到SHD-L2-GIBF 在此條件下比其他算法定位誤差高，即算法的準確度較高。在陣元位置誤差存在下，SHD-L2-GIBF 的誤差均值及方差最小，即算法的魯棒性較好；并且相較最優(yōu)化迭代一范數(shù)方法，此算法步驟計算復雜度較低，物理意義明確，可以量化地衡量聲場，為后續(xù)實現(xiàn)實時處理研究提供較大可能性。

2.4 半消聲室實驗結果

四陣元單球陣的實物圖見圖4(a)，由于實際實驗中硬件設備只有一個四陣元陣，無法實現(xiàn)直接測量32 陣元數(shù)據(jù)，因此實驗時分別對8 個位置測8 次四陣元數(shù)據(jù)，并在每次測試時，于建立的坐標系原點處放置參考傳聲器，通過互相關確定時延從而進行32 陣元采集數(shù)據(jù)的時間同步，得到完整實驗數(shù)據(jù)。

在中國科學院聲學研究所半消聲室進行實驗，實驗搭建場景在圖4(b)中展示。用兩個索尼SRS-XB01 音箱作為聲源，同時發(fā)出等能量的相同白噪聲并用四陣元子陣進行記錄，每次子陣陣列的位置示意圖如圖2(a)中藍色星號所示。兩個音箱分別位于(θ1=?135°, φ1= 75°)和(θ2= 60°,φ2= 120°)處。在300 Hz 處各個算法聲成像結果見圖5。

圖4 子陣及實驗場景實物圖Fig.4 The actual picture of sub-array and experimental scene

圖5 300 Hz 處相干聲源聲成像半消聲室實驗結果歸一化譜圖Fig.5 Normalized spectra of experimental results in a semi-anechoic chamber of coherent sound source acoustic imaging at 300 Hz

根據(jù)聲成像譜圖計算其中3 種可以定位算法在半消聲室實驗中的誤差結果，見表5。

表5 在半消聲室，300 Hz 處兩相干白噪聲聲源3 種算法的定位誤差Table 5 Localization errors of three algorithms for coherent white noise sources at 300 Hz in the semi-0033 cxanechoic room

與2.1 節(jié)仿真結果對比，可見實驗結果與仿真結果類似，說明陣型和算法具有較高的可行性，且SHD-L2-GIBF 算法在此陣列基礎上準確性和魯棒性較好。

3 結論

本文提出了一種基于球諧分解的二范數(shù)廣義逆波束形成(SHD-L2-GIBF)算法，通過在分布式開放球形陣列陣元設置下仿真多種聲成像算法，對比分析可得到此算法可以很好地適配此陣型，在低頻相干聲源的聲成像中提供較好的空間定位能力和算法穩(wěn)定性。在半消聲室中搭建的實驗進一步證明了此算法的可行性。