孫思雨,張海劍,孫 洪
(武漢大學 電子信息學院,湖北 武漢 430072)
聲矢量傳感器(Acoustic Vector Sensor,AVS)是一種由2~3個正交的方向傳感器和一個全向麥克風組成的陣列,方向傳感器用來測量沿正交軸的粒子速度,全向麥克風用來測量該點的聲壓[1-3]。AVS中的各傳感器元件十分靠近,故可假定它們被放置在空間中的同一點上,從而使得AVS具有頻率無關的特性。利用AVS所提供的豐富聲學信息以及其頻率無關的特性,AVS在陣列信號處理方面優于傳統的傳感器陣列。此外,由于具有體積小、質量輕等特點,AVS在聲學研究領域也受到不斷關注,并且在一些新興行業中得到實際應用,如目標跟蹤[4-6]、水聲通信[7]和聲源定位[8]等。
波達方向(Direction of Arrival,DOA)估計是一個備受關注的研究熱點[9-13],隨著AVS在陣列信號處理領域的興起,一些基于AVS的DOA估計方法也紛沓而至。一種常用的方法是將基于傳統壓力傳感器陣列的DOA估計方法擴展到AVS陣列,如波束形成[14-18]、多重信號分類(Multiple Signal Classification,MUSIC)[19]、旋轉不變參數估計(Estimation of Signal Parameters via Rotational Invariance Technique,ESPRIT)[20]以及基于稀疏表示的方法[21]。此外,利用AVS特有的陣列結構,過去三十年中許多學者提出了適用于AVS的聲源測向算法。Nehorai等[22]提出了一種基于強度矢量和速度協方差矩陣的方法,用單個AVS處理DOA估計問題,這是AVS應用在DOA估計中的開創性工作。此外,Levin等[23]還在其中推導了均方角誤差的克拉美羅下界用來評估DOA估計的性能。文獻[24]中提出了基于單個AVS的最大導向響應能量估計器,根據最大功率方向來估計DOA,該方法是基于強度矢量和速度協方差矩陣的擴展,能夠利用合適的參數來接近克拉美勞下界。然而,上述工作主要是為單源測向而設計的,不適用于多源存在的欠定情況。
上述策略大多是在時域內處理觀測信號,而不是充分利用時頻特征,僅利用信號的時域統計信息,難以處理源數超過麥克風數目的欠定情況。隨著時頻分析在信號處理領域的興起[25-26],時頻稀疏性理論逐漸應用于基于AVS的多源DOA估計中。通過檢測單個源占主導的時頻點,并利用這些時頻點進行DOA估計,多源DOA估計任務可以轉換為多個單源DOA估計任務[8]。
準確估計DOA的關鍵在于正確檢測出具有特定特性的有利于測向的時頻點。文獻[27]中提出了一種基于傳感器間數據比(Inner-Sensor Data Ratio,ISDR)的方法,通過提取局域信噪比(Signal to Noise Ratio,SNR)較高的時頻點來實現多源DOA估計[28]。基于ISDR的方法在自由場假設的前提下取得了較好DOA估計性能,但隨著混響的增加,反射信號會破壞ISDR,從而導致該方法整體性能下降。為提高在混響環境下的DOA估計性能,一般選擇受混響干擾較少的時頻點從而實現精確DOA估計。文獻[29]中提出了一種基于低混響單源點(Low-Reverberant-Single-Source Point,LRSSP)的方法,利用單個AVS在混響環境中實現多源DOA估計,其中LRSSP指代時頻域中受混響干擾較小的單源時頻點。針對AVS特有的頻率無關特性,該方法利用每個LRSSP的實部和虛部具有相似方向向量的特點,用來在時頻域中檢測LRSSP。盡管基于LRSSP方法具有較好性能且易于實現,但隨著混響增加,該方法不可避免地會檢測到一些包含錯誤方向信息的非LRSSP從而導致性能下降。
為解決上述估計問題,進一步提高強混響情況下DOA估計的準確性和魯棒性,本文提出了一種有效的方法用來檢測直達聲單源點(Direct-Path-Single-Source Point,DP-SSP),指代由單個實際源信號而非反射源所占主導的時頻點。DP-SSP中包含了相關源的正確方向信息,能防止估計的方向與真實DOA產生較大偏差。本文受文獻[30]直接路徑優勢(Direct Path Dominance,DPD)測試的啟發,提出的DP-SSP檢測技術首先使用頻率平滑來減輕反射信號的影響,然后利用子空間投影技術來進一步提取DP-SSP并消除干擾項。檢測出的DP-SSP隨后被用來進行源數目估計和聚類。最后將MUSIC方法用于每一類DP-SSP點簇中用來計算該點簇對應源的DOA。下文將詳細介紹本文所提出的基于頻率平滑和子空間分解的DP-SSP檢測方法。
假設一個密閉混響室內中存在M個聲源和一個由3個方向傳感器和一個全向麥克風組成的AVS陣列,該AVS陣列接收到的混合信號經過短時傅里葉變換(Short-Time Fourier Transform,STFT),得到時頻域信號可表示為:
(1)
式中:t和f表示分別表示時間軸和頻率軸索引,X(t,f)=[Xp(t,f),Xvx(t,f),Xvy(t,f),Xvz(t,f)]T中包含了AVS四個通道信號的STFT系數,Xp(t,f)、Xvx(t,f)、Xvy(t,f)、Xvz(t,f)分別表示全向麥克風和3個方向麥克風的輸出,Hq(f)表示第q個源到AVS之間的聲傳輸函數,Sq(t,f)表示第q個源的STFT系數,N(t,f)表示噪聲的STFT系數。
由于混響環境下的觀測信號一般由直接路徑聲和反射聲組成,因此第q個源到AVS之間的聲傳輸函數可以寫為:
(2)


圖1 信號反射模型Fig.1 Signal reflection model
這里Iq表示強反射信號等效產生的虛擬鏡像源的個數,ai(f)表示第i個源到AVS之間的聲傳遞函數,Ψi=(θi,φi)表示第i個源的DOA,則基于AVS的陣列流形第i個源的導向矢量可以表示為[24]:
d(Ψi)=[1,cosφicosθi,sinφicosθi,sinθi]T,
(3)
以AVS所處位置為笛卡爾坐標系原點,φi∈[0°,360°)和θi∈(-90°,90°]分別表示第i個源到原點的方位角和仰角。
基于上述分析,本文的鏡像模型可以等效為一個不含混響的遠場模型:
(4)
式中:J表示實際源和虛擬源總數,其中包含Q個實際源和J-Q個由強反射產生的虛擬鏡像源,E(t,f)表示建模誤差,包含弱混響或者低能量的干擾源點。因此,在混響環境下的信號模型可以近似轉化為自由空間中存在J個源的遠場模型。在低混響情況下,J的值與Q接近。
根據式(4)中描述的信號模型,可以將時頻點分為DP-SSP、強反射單源點(Strong Reflection Single Source Point,SR-SSP)以及多源時頻點(Multi-Source Point,MSP)。DP-SSP是一種能量由一個實際源的直接路徑信號占主導,受其他信號源和混響影響較少的時頻點。SR-SSP表示由強早期反射信號等效的鏡像源占主導的時頻點,強早期反射信號通常指第一次反射的信號。MSP表示由多個能量值較大的源占主導的時頻點,其中包含實際源和虛擬鏡像源。由于SR-SSP和MSP會對最終的DOA估計結果產生不利影響,算法設計時需要考慮去除這類離群值。
從理論上來說,少量的DP-SSP已經足以獲得準確的DOA估計。然而DP-SSP的檢測通常會伴隨著一系列的干擾時頻點,這些時頻點在檢測過程中不可避免地會被錯誤識別成DP-SSP。因此,為了精確估計DOA,顯然希望正確的DP-SSP點在所有檢測出的時頻點中占比越高越好。文獻[29]提出了基于LRSSP的方法,致力于利用AVS接收到的信號存在的獨特數學特性來設計算法提取DP-SSP,從而獲得準確的源測向結果。然而隨著混響增加,大量滿足LRSSP檢測規則的異常值被錯誤檢測為DP-SSP,最終導致真實的DP-SSP在所有檢測出的時頻點中占比下降,進而對最終DOA估計造成影響。
為了分析基于LRSSP[29]的方法在高混響條件下識別DP-SSP的性能,給出了每個檢測到的LRSSP點的強度向量分布,并計算其中真實DP-SSP的占比。檢測到的LRSSP的強度向量計算方法如下[29]:
(5)
(6)
式中:{·}*表示共軛操作,uq表示第q個源的單位方向向量,當強度向量η(t,f)近似與一個源的單位方向向量平行時,則其對應的時頻點(t,f)被認為是該源的一個DP-SSP點,其判別規則可以表示為:
(7)

(8)
式中:Card(·)表示一個集合中元素的數量,集合Ωd中包含了所有符合檢測規則的時頻點,Ωd中所有滿足式(7)中的時頻點包含在集合ΩrD中。在本文中,實際DP-SSP占比P被用來作為評估DP-SSP檢測性能的一個標準。
本文實現一個簡單的仿真來輔助說明上述分析,對比結果展示在圖2中,不同顏色的圓圈表示對應于不同源的真實DP-SSP。仿真房間大小為8 m×6 m×4 m,AVS被放置在房間中心,3個聲源來自于TIMIT數據集[32],均被放置于離聲矢量傳感器1.5 m處,3個源的DOA分別為(50°,0°)、(110°,15°)和(170°,30°)。

(a)基于LRSSP的方法(RT60=0.3 s)

(b)基于LRSSP的方法(RT60=0.6 s)圖2 檢測到的時頻點與真實DP-SSP點對應的強度矢量分布圖(SNR=20 dB)。Fig.2 The intensity vector distribution of the detected time-frequency points and real DP-SSP (SNR=20 dB)
從圖2中可以明顯看到,檢測到的LRSSP中僅有一小部分是對DOA估計有利的真實DP-SSP點。通過對比圖2(a)和圖2(b)可以發現,隨著混響的增加,強度向量分布的方向特征變得不明顯,且實際DP-SSP的比例逐步下降。盡管基于LRSSP[29]的方法會對檢測到的時頻點做進一步的離群值去除使強度矢量分布具有更加明顯的方向性,但是這僅適用于低混響的情況。隨著混響的增加,大量出現的離群值會使得每個源對應的強度矢量簇中心與真實的方向向量發生偏移,進而在進行離群值去除時導致大量有利于DOA估計的時頻點被去除,從而影響最終DOA估計的精度。為解決上述實際DP-SSP占比較低的問題,本文提出了一種有效的DP-SSP檢測方法,通過準確檢測DP-SSP并消除離群值,來提升有利DP-SSP的整體占比。
基于上述分析可以得到,在高混響環境下,DP-SSP的檢測對DOA估計來說是至關重要的,因為DP-SSP中包含了對應源的關鍵方向信息。隨著混響的增強,DP-SSP的檢測變得越來越有挑戰。本文提出了一種有效的兩階段策略來檢測有用的DP-SSP并去除干擾時頻點,在第一階段應用頻率平滑技術來盡可能去除潛在的早期反射強源點SR-SSP,在第二階段應用子空間投影進一步有效去除多源點MSP并提取DP-SSP。
2.2.1 采用頻率平滑消除SR-SSP
早期反射信號被認為是對直達聲信號做延遲和衰減。隨著混響增強,早期反射信號的能量逐漸增加變得不可忽略,因此第一次反射的信號通常被建模為一個與實際源信號高度相關的虛擬鏡像源發出的信號。SR-SSP指代那些由虛擬源占主導的時頻點。為了減少SR-SSP對最終DOA估計的影響,首先對所有的時頻點做一次頻率平滑以消除潛在的SR-SSP點,并粗略提取不受混響干擾的DP-SSP[33]。
受DPD測試[30]的啟發,本文對分別對語譜圖在時間和頻率軸上平滑取平均,構建局部時頻相關矩陣:

(9)
式中:Jt和Jf分別表示滑動窗口的尺寸大小。鑒于AVS的頻率無關特性,導向矢量在頻率平滑的過程中不會受到影響,因此上述時頻相關矩陣可表示為:
(10)


(11)

(12)
式中:σi(t,f),i=1,2,…,J表示矩陣σi(t,f),i=1,2,…,J的一個特征值。如果第i個源在該時頻區域內不活躍,則其對應的特征值將接近于0。將式(12)帶入式(10)得到:
σJ(t,f)d(ΨJ)Td(ΨJ)。
(13)

(14)
因此,潛在DP-SSP的檢測規則可以寫為:

(15)
式中:σmax和σsubmax分別表示最大和第二大的特征值。
通過頻率平滑和特征值分解操作去除了潛在的SR-SSP,減輕了混響的影響。由于上述操作是一個“區域級”的操作,會不可避免地檢測到許多非DP-SSP點,比如多源點MSP,因此本文將第一階段檢測得到的點重新定義為低混響自項源點(Low Reverberant Auto Source Points,LR-ASP)。所有滿足式(15)中檢測規則的時頻點被包含在集合ΩA中。為了進一步精確提取DP-SSP,將在下一小節對第一階段檢測到的LR-ASP進行“點級”的處理。
2.2.2 采用子控件分解消除MSP
由于集合ΘA中的LR-ASP主要包含DP-SSP和MSP,因此需要對LR-ASP進行提純,在“點”級別上提取DP-SSP點。首先利用信號時頻稀疏性,對每個LR-ASP做如下假設:在每個LR-ASP出現的強能量源點不超過2個。這個假設在大多數情況下是可行的,因為涉及3個或更多高能量源的時頻點數量太少可忽略不計。
因此,每個LR-ASP可以寫作:
(16)
式中:q1,q2∈{1,2,…,Q}表示出現在該時頻點上的2個源的序號。當2個源之間能量值比值大于一個遠大于1的閾值時,能量值高的源被認為主導該時頻點,可以被看作是對應源的一個DP-SSP點,該閾值一般為經驗選取。計算并歸一化每個時頻點對應的空間向量可得:
Cq1(t,f)d(Ψq1)+Cq2(t,f)d(Ψq2),
(17)


(18)
然而,當LR-ASP是一個由2個源占主導的MSP,上述條件將不會被滿足。接著對所有的歸一化空間向量進行k-means聚類,得到N簇向量,每個簇心對應的向量構成了一個過完備導向矢量字典。由于歸一化空間向量是一個方向向量,對應的是與單個源相關的DP-SSP或者與多個源相關的MSP,因此將N設置為一個大于Q的值。每個簇心對應的方向向量計算為:
(19)

(20)

根據前述假設,每個LR-ASP上活躍的源信號不超過2個,則在時頻點(t,f)上最可能出現的2個源的STFT系數可以計算為:
(21)
式中:[·]?表示Moore-Penrose偽逆操作符,aq1、aq2表示2個源對應的導向矢量。在字典ΘA中搜索,找到最優的2個導向矢量aq1、aq2,滿足下式:
(22)
(23)
A2=[ar1ar2],
(24)
式中:P表示矩陣A2噪聲子空間上的正交投影矩陣,ar1和ar2分別表示從字典ΘA中隨機挑選的2個列向量,I表示單位矩陣,(·)-1表示求逆操作。
2個源的STFT系數可由式(21)計算得到。通過計算2個源之間的能量比,并設定閾值,即可挑選出一組DP-SSP點:

(25)
式中:λ2為一個經驗閾值。所有滿足上述不等式的LR-ASP點都被認為是DP-SSP點,包含在集合ΩD中。通過對集合ΩA中的所有LR-ASP點進行子空間分解,可以有效去除2個強能量源占主導的多源點MSP。
為了驗證所提的DP-SSP檢測方法的有效性,首先在對基于LRSSP[29]的方法應用頻率平滑,去除潛在的SR-SSP點。結合LRSSP[29]檢測和SR-SSP去除后得到的時頻點強度矢量分布如圖3所示。雖然最后檢測到的DP-SSP分布的方向性仍不明顯,但與原始基于LRSSP[29]檢測的方法相比,真實DP-SSP所占比例有明顯提高。本文所提的兩階段DP-SSP檢測方法檢測到的DP-SSP強度矢量分布如圖3(b)所示,其分布具有明顯的方向性。雖仍然存在一些方向向量偏離真實DOA的較為分散的離群點,但總體強度矢量的分布非常接近真實源方向。此外,采用所提出的DP-SSP檢測方法,相比基于LRSSP[29]的方法,真實DP-SSP占比大大提高,有力地說明了該方法在檢測DP-SSP上的優越性。

(a)LRSSP檢測+SR-SSP去除 (RT60=0.6 s)

(b)本文所提DP-SSP檢測算法 (RT60=0.6 s)圖3 2種檢測規則檢測到的時頻點與真實DP-SSP點對應的強度矢量分布圖(SNR=20 dB)Fig.3 The intensity vector distribution of the detected time-frequency points and real DP-SSP using different detection rules (SNR=20 dB)
隨后利用平滑直方圖[29]方法對檢測到到DP-SSP進行源數估計,然后根據估計出的源數目將所有檢測到的DP-SSP聚類為多組。最后對每組DP-SSP利用MUSIC方法進行DOA估計,整合每組DOA估計結果得到最終多源DOA估計結果。
給出了一系列實驗結果用來評估所提算法在仿真環境中的DP-SSP檢測、源數估計以及DOA估計3個方面的性能,并與基于LRSSP[29]的方法進行了比較。為了突出所提兩階段DP-SSP方法每一步的有效性,本文還將SR-SSP去除與LRSSP[29]檢測結合作為另外一種對比方法,用來評估DP-SSP的檢測性能。
仿真混響室內為一個理想的封閉長方體,尺寸為8 m×6 m×4 m,AVS被放置在房間中心。采用鏡像原理[31]生成混響室內的房間脈沖響應,聲源被放置在以AVS為圓心,半徑1.5 m的圓上。聲源從TIMIT語料庫[32]中挑選,采樣率為16 kHz。麥克風接收到的信號進行STFT變換,FFT點數為1 024,幀間重疊為50%。加性噪聲為高斯白噪聲,其他實驗參數設置如表1所示,表1中的閾值均為經驗選取。

表1 參數設置Tab.1 Parameter settings
與文獻[29]類似,本文采用均方角誤差(Root-Mean-Square Angular Error,RMSAE)作為評估指標來量化所提方法在所有實驗中的整體性能,RMSAE定義為:
(26)
(27)

不同混響環境下利用3種不同的檢測規則檢測到的時頻點中真實DP-SSP占比如圖4所示。

圖4 采用3種方法得到的時頻點中真實DP-SSP所占百分比隨混響時間變化Fig.4 Variation of real DP-SSP percentage with reverberation using three detection rules
從圖3中可以看到,隨著混響的增強,基于3種檢測規則的真實DP-SSP占比呈現下降趨勢。但隨著混響增加,所提的兩階段DP-SSP檢測方法始終保持著明顯優勢。結合SR-SSP去除和LRSSP[29]檢測的方法在DP-SSP檢測上優于原始的LRSSP[29]方法,這說明對去除SR-SSP可以有效輔助DP-SSP的檢測。通過對比也可發現,除了第一階段的SR-SSP去除,第二階段的子空間投影算法也在DP-SSP檢測中展現極大優勢。真實DP-SSP占比的提高使得所提算法能在高混響情況下實現更高精度的源計數和DOA估計。
表2給出了所提DP-SSP檢測方法與基于LRSSP[29]的方法相比,不同混響環境下正確估計源信號數目的概率。從表中可以發現,2種檢測規則在RT60=0.3 s的情況下均能準確估計源數,信噪比對估計結果的影響較小。當RT60<0.45 s,SNR=20 dB時,基于LRSSP[29]的方法基本可以正確估計源數,準確率大于90%,且當RT60<0.45 s時,基于LRSSP[29]的方法估計精度斷崖式下降。此外,從表中可以看到,基于LRSSP[29]的方法性能受SNR影響很大,這說明該方法對噪聲更加敏感。本文所提方法隨著混響增強和信噪比降低,仍然能保持穩定有效的源數估計性能。

表2 在不同信噪比和混響程度下正確估計源數目的概率Tab.2 Percentage of the correct estimation of the number of sources with different SNRs and reverberation
圖5中展示了SNR為20 dB,在RT60=0.6 s和RT60= 0.75 s的仿真混響環境下,不同算法進行DOA估計的可視化結果。3個源信號所在位置分別為(70°,0°)、(110°,15°)、(170°,30°)。然而,當混響時間增加到0.75 s時,基于LRSSP[29]的方法無法準確定位出每個源,如圖5(b)所示,一些源的DOA估計與對應的真實DOA有很大偏差。與基于LRSSP[29]的方法不同,隨著混響增加,所提算法仍能實現高精度的DOA估計,且誤差較小。此外,圖6展示了在多個源間隔較近的情況下DOA的估計結果。3個源的DOA分別為(70°,0°)、(110°,15°)、(150°,25°)。結果表明,基于LRSSP[29]的方法無法處理距離較近的多源DOA估計,而所提方法能實現準確估計,具有明顯優勢。

(a)基于LRSSP[29]的方法(RT60=0.60 s) (b)所提DP-SSP檢測的方法(RT60=0.60 s)

(c)基于LRSSP[29]的方法(RT60=0.75 s) (d)所提DP-SSP檢測的方法(RT60=0.75 s)圖5 仿真環境下2種方法估計DOA的可視化結果圖(SNR=20 dB)Fig.5 DOA estimation results in the simulated environment using different DOA estimation algorithms (SNR=20 dB)

(a)基于LRSSP[29]的方法(RT60=0.60 s) (b)所提DP-SSP檢測的方法(RT60=0.60 s)圖6 源相隔較近時2種方法估計DOA的可視化結果圖(SNR=20 dB)Fig.6 DOA estimation results using different DOA estimation algorithms in the case of closely spaced sources (SNR=20 dB)
表3展示了利用2種不同算法計算的RMSAE隨混響時間變化的情況。由表3可以看出,在RT60<0.6 s時,基于LRSSP[29]的方法在3個源DOA估計上表現較好,但當混響時間超過0.6 s時,該方法的DOA估計精度明顯下降。值得注意的是,盡管所提方法的性能不可避免地隨著混響的增加而下降,但當RT60增加到0.9 s時,在SNR=20 dB的條件下,對比算法已經完全無法進行多源DOA估計,然而所提算法仍然具有一定性能。當信噪比降到15 dB,且RT60=0.9 s時,所提方法性能略有下降。

表3 RMSAE隨混響和SNR變化Tab.3 Variation of RMSAE with the reverberation and SNR
圖7展示了對不同數量的源進行DOA估計后計算得到的RMSAE隨混響變化的趨勢。結果表明,隨著源數量的增加,本文所提方法性能不可避免地下降,特別是存在4個源的情況下,性能下降明顯,這是因為DP-SSP檢測的準確性隨著源數量的增加而下降。所提算法在2個源的情況下,在不同的環境中都能展現出較好的性能,即使當RT60增加到0.9 s時,仍然能獲得較高精度。

圖7 不同源數情況下使用所提方法計算的RMSAE隨混響時間變換圖(SNR=20 dB)Fig.7 Variation of RMSAE with reverberation time in case of different number of sources by the proposed algorithm (SNR=20 dB)
本文針對傳統算法利用單個AVS在高混響情況下多源DOA估計性能差的問題,提出了一種新的多源DOA估計算法。該算法的關鍵在于檢測DP-SSP,這些DP-SSP為DOA的估計提供了有利信息。本文所提出的DP-SSP檢測算法有兩方面的貢獻。首先,針對單個AVS設計了一種基于頻率平滑的SR-SSP去除方法,有效地消除了混響帶來的不利影響,為準確檢測DP-SSP提供幫助。其次,采用子空間投影的算法進一步細化DP-SSP的檢測,消除大部分MSP,提取正確的DP-SSP。實驗結果表明,該算法對DP-SSP的高精度檢測使得本文算法在源數目估計和DOA估計方面取得較好性能。