中國人標準頭模BHead210的頭相關傳遞函數(shù)數(shù)值計算

唐玲，付中華

?

中國人標準頭模BHead210的頭相關傳遞函數(shù)數(shù)值計算

唐玲1，付中華2

(西北工業(yè)大學計算機學院，陜西西安 710129)

針對基于中國人平均頭部尺寸的BHead210標準頭模，通過結構光掃描獲取三維數(shù)據(jù)并利用邊界元方法計算其頭相關傳遞函數(shù)。將剛球模型的計算結果與解析結果進行對比，驗證了計算方法的正確性，并將BHead210的計算結果與測量結果對比以觀察誤差。給出了由邊界元計算數(shù)據(jù)導出的BHead210雙耳時間差、雙耳聲級差曲線與不同聲源距離某方位的頭相關傳遞函數(shù)幅度譜。最后對比和分析了BHead210頭模數(shù)據(jù)與CIPIC公布的KEMAR頭模數(shù)據(jù)的差異。

頭相關傳遞函數(shù)；中國人標準頭模BHead210；數(shù)值計算

0 引言

頭相關傳遞函數(shù)(Head Related Transfer Function, HRTF)定義為自由聲場下從聲源到雙耳的頻域聲學傳輸函數(shù)，它表達了生理結構對聲波的綜合濾波效果[1]。在虛擬聲的研究領域中，尤其在雙耳空間聽覺研究方面，HRTF起著重要的作用。

獲取HRTF的方法有實驗測量和理論計算兩種。實驗測量可在消聲室中進行，采用的是數(shù)字測量技術[2]。實驗測量對實驗環(huán)境和設備有很高的要求，需要大量的校準和均衡處理，是一項非常耗時的工作。理論計算則是通過求解頭部、軀干、耳廓等組成的聲學系統(tǒng)對聲波的散射問題來獲取HRTF[1]。近20多年來，許多研究工作者開展了理論計算HRTF的工作。1998年，Kahana等人采用邊界元方法對KEMAR人工頭HRTF的計算頻率上限達到了6kHz，并且證明了與實驗結果相符合[3]。2003年，TIAN Xiao等用時域有限差分法結合完全匹配層方法計算了球模型和KEMAR人工頭的HRTF[4]。2004年，F(xiàn)els等人為了研究HRTF隨兒童生長的變化，通過數(shù)字照相獲取兒童(4~6歲)頭部三維信息并用邊界元方法計算了兒童的HRTF[5]。2007年，Wolfgang等人應用由傳統(tǒng)邊界元發(fā)展而來的多級邊界元方法計算HRTF，大大提高了理論計算的速度[6]。

目前，國外已經(jīng)公開部分的人工頭及真人的HRTF數(shù)據(jù)庫，如MIT的Gardner等于1994年對KEMAR人工頭進行測量獲得的HRTF數(shù)據(jù)庫[7]，加州大學Davis分校CIPIC實驗室的Algazi 和Duda等于2000年通過測量獲得的真人與KEMAR人工頭的HRTF數(shù)據(jù)庫[8]。國內(nèi)華南理工大學謝菠蓀教授等于2005年底也通過對52名受試者的測量獲得了相應的中國人頭的HRTF數(shù)據(jù)庫[1]。

HRTF是具有個性化特性的物理量，它與人體頭部、耳廓、軀干等生理結構有著密切的聯(lián)系，國外的數(shù)據(jù)庫是基于西方人的生理特征而建立的，然而中國人在生理外形和尺寸上與西方人有一定差別，所以國外數(shù)據(jù)庫不一定完全適用于中國人[1]。中國傳媒大學依據(jù)中國成年人頭面部尺寸的國家標準設計與制作了仿真頭模BHead210[9]，獲取BHead210的HRTF工作也剛剛展開。本文基于邊界元法，通過對BHead210掃描獲取三維信息計算其HRTF。

1 邊界元法計算HRTF

1.1 原理簡介

邊界元法(Boundary Element Method, BEM)是以邊界積分方程為控制方程，對邊界進行分單元離散，再將所有單元上的控制積分方程整合為代數(shù)積分方程組進行求解。

邊界元分為直接邊界元與間接邊界元，直接邊界元的計算區(qū)域只能在邊界腔體的內(nèi)部或外部，而間接邊界元的計算區(qū)域則是整個空間[10]。邊界元方法的一個關鍵問題是，對于某些頻率會產(chǎn)生奇異值，這是因為線性方程組的求解矩陣不滿足非奇異條件使得在該頻率不存在唯一的解。直接邊界元方法可以應用亥姆霍茲積分公式方程組(Combined Helmholtz Integral Equation Formulation, CHIEF)[11]方法解決解的非唯一性問題，盡管如此，間接邊界元方法比結合了CHIEF的直接邊界元法能更好地處理奇異問題[12]。

若分別計算不同距離、不同方位聲源的HRTF將耗費許多時間，使用聲學互易原理可以大大減少計算量。聲學互易原理指在線性聲學范圍內(nèi)，從發(fā)射點到接收點之間的聲學系統(tǒng)是一個互易系統(tǒng)[13]。

本文采用LMS公司的Virtual.Lab Acoustic軟件(簡稱VL)[14]中的邊界元模塊中的間接邊界元方法來計算HRTF數(shù)據(jù)。

1.2 中國人工頭模型及掃描

中國人標準頭模BHead210是按照中國成年人頭面部尺寸國家標準GB/T 2428—1998設計與制作的仿真頭模(如圖1所示)，其頭形參數(shù)、結構均比較符合中國人的特點[9]。在該頭模上安裝的耳廓模型，是在對中國成年男女耳廓形態(tài)進行了群體研究的基礎上，選取最具有代表性的耳廓外形結構作為參考而設計制作的“平均耳”模型(或標準耳廓)[9]。該頭模在耳道入口處放置了兩個微型麥克風，因此可以通過實驗測量獲取HRTF數(shù)據(jù)并與數(shù)值計算結果進行比較分析。

圖1 BHead210標準人工頭模

在計算HRTF之前，首先要獲取人工頭模的三維掃描結果，并合理劃分表面網(wǎng)格。下面介紹建立頭模網(wǎng)格模型的方法，主要分為三維掃描與網(wǎng)格處理兩步。

首先用掃描方法獲取頭模的三維結構信息。通常可以采用激光三維掃描儀或結構光三維掃描儀進行掃描。本文使用DAVID結構光三維掃描儀[15]來獲取頭模信息。結構光三維掃描儀相對于激光三維掃描儀有設備便攜、成本低等優(yōu)點[15]。結構光三維掃描儀一次只能掃描到光直射的范圍，因而需要對頭模的各個面分開進行掃描。掃描完成后，將所有面融合成一個完整的人工頭模型。融合的結果可能會在交接處不光滑，可以利用三維軟件將之平滑。掃描得到的結果是以obj文件格式保存的，在該文件里包含了組成頭模的三角形面網(wǎng)格信息，約有63萬個三角形網(wǎng)格，平均邊長為0.8 mm。

在邊界元計算中，分析頻率上限的條件是：模型最大的網(wǎng)格邊長應在所要計算的最短波長的1/4~1/6之間[16]。若計算的上限頻率為10 kHz，則網(wǎng)格的最大邊長約為6.7 mm。通過掃描得到的網(wǎng)格模型并不理想，不僅網(wǎng)格數(shù)量過多增大了計算量，而且由于融合過程中分辨率的不同導致網(wǎng)格大小可能不均勻，在光滑的部分網(wǎng)格邊長過大，在精細部分邊長又過小，如頭頂?shù)木W(wǎng)格邊長相對于耳朵要大得多，因此需要對網(wǎng)格模型進行處理。首先將掃描得到的網(wǎng)格模型進行從小到大的合并，減少網(wǎng)格數(shù)量。本文采用了rapidform軟件[17]的decimate功能對掃描結果進行網(wǎng)格合并。網(wǎng)格合并得到的新網(wǎng)格模型保留了原網(wǎng)格模型的特點：平滑處的網(wǎng)格邊長較大，非平滑處的網(wǎng)格邊長很小，其中耳朵處的結構最為復雜，所以該處的網(wǎng)格最小(如圖2所示)。網(wǎng)格合并時保證耳朵處網(wǎng)格邊長為4 mm，其他部位的網(wǎng)格均大于4 mm。接著需要對平滑處過大的網(wǎng)格再細分，使其滿足計算條件。細分方法是，將邊長大于6.5 mm的三角形網(wǎng)格沿該邊長對半分成兩個三角形，直到所有邊長均小于6.5 mm為止(結果如圖3所示)。最后，將新的網(wǎng)格模型在rapidform中再進行網(wǎng)格的重新劃分，本文利用了rapidform中的remesh功能[17]。經(jīng)過上面的重新整理，新的網(wǎng)格模型邊長都在4~6.5 mm之間，滿足HRTF的計算要求，總單元格數(shù)為2.5萬(如圖4所示)。網(wǎng)格的處理過程中，軟件并未改變頭模形狀，所以處理后的結果將不影響計算結果。

圖2 處理前耳朵處網(wǎng)格與頭頂處網(wǎng)格對比

圖3 處理后耳朵處網(wǎng)格與頭頂處網(wǎng)格對比

圖4 人工頭網(wǎng)格模型

2 計算結果與分析

剛球模型是分析HRTF數(shù)據(jù)的最簡單模型，為了驗證邊界元計算方法的正確性，將剛球模型的邊界元法計算結果與解析結果進行了對比并分析誤差。隨后給出了中國人標準頭模BHead210的HRTF相關數(shù)據(jù)，包括邊界元法計算結果與實際測量結果的部分方位對比圖、由計算得到的HRTF數(shù)據(jù)導出的雙耳時間差與雙耳聲級差曲線、某方位不同聲源距離的HRTF幅度譜等。實驗測量目前仍然是獲取HRTF數(shù)據(jù)最為準確的方法，因而有必要將計算與測量結果進行對比。雙耳時間差()、雙耳聲級差()是聲源的定位重要因素，可通過HRTF數(shù)據(jù)計算得到。此外，理論計算相對于實驗測量，還有容易獲取近場HRTF的特點，可方便地計算不同聲源距離的HRTF以分析HRTF與聲源距離的關系。最后將BHead210頭模的計算結果與CIPIC公布的KEMAR頭模數(shù)據(jù)進行了對比，觀察BHead210數(shù)據(jù)與KEMAR數(shù)據(jù)的差異。

2.1 剛球模型HRTF

剛球模型是略去了耳廓、軀干等作用，把人頭簡化為半徑的剛性球體，瑞利是最早采用剛球模型分析雙耳聲壓的。Duda等利用球體對點聲源產(chǎn)生的球面波的散射公式，得到了球體模型的HRTF解析公式[18]：

取剛球體半徑為87.5 mm，聲源位置離球中心為1.5 m，分別用邊界元方法與式(3)計算HRTF，計算頻率范圍為100~10000 Hz，步長為70 Hz，將得到的結果進行對比(如圖5所示)。圖5中角度為球中心到球面觀測點(耳朵位置)的直線與聲波入射方向之間的夾角，如0°表示聲源與觀測耳同側，90°表示正前方，180°表示聲源與觀測耳異側。

在低頻處，邊界元計算結果與理論結果相差甚小，隨著頻率增高，誤差會相對大一些。此外，聲源與耳朵同側時(如0°曲線)，邊界元計算結果誤差比較小，而當聲源與耳朵異側時(如150°、180°曲線)，計算誤差會相對增大。得到的實驗數(shù)據(jù)中，最大絕對誤差不超過1.2 dB，所以我們采用的邊界元法計算剛球模型HRTF誤差是可以接受的。

圖5 球體HRTF

2.2 中國人工頭HRTF

本節(jié)給出了HRTF計算結果與測量結果的對比圖(見圖6)，同時為了觀察HRTF數(shù)據(jù)的基本特性，通過邊界元計算得到的HRTF數(shù)據(jù)導出了雙耳時間差()、雙耳聲級差()曲線(如圖7、8所示)，并計算了不同聲源距離下某方位的HRTF幅度譜(圖9、10)。

圖6是水平面?zhèn)€別方位上BHead210的HRTF數(shù)據(jù)邊界元計算與實驗測量結果對比圖，其中實驗測量結果來自中國傳媒大學的傳播聲學研究所[9]。將頭模中心設置在連接兩耳線段的中心點，點聲源放置于離頭模中心距離1 m處，觀察點為耳道入口處。在邊界元法計算中，假設頭模為理想的剛體模型，計算頻率范圍為0~12 kHz，步長為70 Hz。圖中分別給出了仰角0°(即水平面)，方位角0°(正前方)、90°(正右方)、180°(正后方)、270°(正左方)的左耳HRTF幅度譜。

從圖6可以看出數(shù)值計算方法和實驗測量結果比較接近，曲線走勢基本一致。進一步觀察數(shù)據(jù)結果時發(fā)現(xiàn)，在高頻部分，絕對誤差會有所增大，而且部分方位的計算與實驗測量結果還存在較大的偏差(如觀測耳的異側方向)。耳廓谷的位置和幅度值在某些方位有明顯差異，正前方差異較大，正后方與觀測耳的同側方向的差異在1 kHz以內(nèi)。數(shù)據(jù)上造成的不一致性可能來自于頭模的剛性假設與掃描誤差，在主觀聽覺上也需要進一步地探討。

為了進一步驗證計算結果，我們將計算得到的HRTF數(shù)據(jù)導出了、曲線。是聲源方向定位的一個主要因素，在頻率小于約1.5 kHz的低頻段，它對方向定位起著重要的作用。在頻率約為1.5~4 kHz的情況下，、對方向定位共同起作用。在頻率大于4~5 kHz的高頻段，是方向定位的主要因素。

圖6 BHead210計算及測量結果對比圖

圖7為由計算得到的HRTF導出的曲線，本文采用的是頻域的相關法計算，得到的結果是全頻帶范圍的“平均”。如圖7所示，在方位角0°與180°處為0，隨著聲源接近側向，增加，在側向附近達到最大。

圖8為不同頻率下的曲線，在低頻處很小，且隨方位角的變化也很小。隨著頻率增大，曲線也變得復雜，如1.4 kHz、2.8 kHz的曲線最大值并不在90°方向，這是由于聲波通過頭的前、后、上方繞射并相干疊加到達異側耳而形成的聲壓亮點。

圖7 ITD曲線

圖8 ILD曲線

為了觀察聲源距離與HRTF數(shù)據(jù)的關系，圖9、10給出了聲源處于不同距離下，水平面方位角分別為90°、270°的左耳HRTF幅度譜，計算條件與前面所提的一樣。實驗測量獲取近場HRTF數(shù)據(jù)是有一定困難的，理論計算則可以容易地獲取不同聲源距離的HRTF數(shù)據(jù)。從圖中可看出，當聲源距離小于0.5 m時，幅度譜的變化隨距離的變化比較大。隨著聲源距離的增加，幅度譜的變化也逐漸變小。因而，在聲源距離大于1~1.2 m時，HRTF的測量與分析可作為遠場處理，近似于距離無關。

2.3 中國人工頭與KEMAR人工頭對比

HRTF數(shù)據(jù)是與人體頭部、軀干、耳廓等外形尺寸和結構有關的具有個性化特征的物理量。國外公開的HRTF數(shù)據(jù)是以外國人特征為基礎獲得的。BHead210頭模是依據(jù)中國人的平均標準特征而設計，我們將計算得到的BHead210的HRTF數(shù)據(jù)與CIPIC公布的KEMAR人工頭模HRTF數(shù)據(jù)進行對比，以觀察兩者的差異。

圖9 90°左耳不同距離HRTF幅度譜

圖10 270°左耳不同距離HRTF幅度譜

圖11為聲源距離在1m處產(chǎn)生的水平面0°(正前方)、90°(正右方)、180°、270°的中國人工頭BHead210與KEMAR人工頭的右耳HRTF數(shù)據(jù)結果圖。BHead210的數(shù)據(jù)為邊界元計算數(shù)據(jù)，計算條件與前面所提一樣。KEMAR人工頭配有大小號耳廓，數(shù)據(jù)為測量數(shù)據(jù)，觀測點在耳道入口處[8]。KEMAR人工頭的說明書建議小號耳廓適用于歐美女性和日本人，而大號耳廓適用于歐美的男性。

觀察對比圖，KEMAR的大小耳廓的HRTF值在低頻段是相近的，只有在高頻處才出現(xiàn)明顯差異，耳廓谷的方位和大小也有所差異。相比之下，KEMAR與BHead210的結果在整個頻段都有一定差異。在低頻段，KEMARA與BHead210在正前方0°位置差異較大。在高頻段，BHead210與KEMAR差異也比較明顯，尤其是異側耳處的耳廓谷方向。這些差異在主觀聽覺上的影響還有待進一步研究。

圖11 中國人工頭與KEMAR的HRTF對比圖

3 結 論

HRTF包含了聲源的空間定位信息，它是一個具有個性化特性的物理量，通過實驗獲取HRTF是一項復雜的工作，理論計算可使工作得到簡化。獲取中國人標準頭模BHead210的HRTF數(shù)據(jù)，可為分析中國人虛擬聽覺特性提供一些數(shù)據(jù)基礎。本文介紹了獲取中國人標準頭模BHead210的HRTF數(shù)據(jù)的理論計算方法，首先通過結構光掃描和網(wǎng)格整理獲得了BHead210的三維網(wǎng)格結構模型，然后采用邊界元方法結合互易原理計算出不同方位的HRTF數(shù)據(jù)。對比了剛球的解析結果和邊界元計算結果，驗證了邊界元方法的正確性，然后將中國人標準頭模BHead210的計算結果與實際測量結果進行誤差分析，并給出了由計算的HRTF數(shù)據(jù)導出的、曲線以及不同聲源距離的HRTF幅度譜以分析HRTF特性，最后還將計算得到的BHead210人工頭的HRTF與CIPIC公布的KEMAR人工頭的HRTF數(shù)據(jù)進行對比。

致謝：作者在此衷心地感謝來自中國傳媒大學的傳播聲學研究所的齊娜老師，她為本文提供了中國標準頭模BHead210的實驗測量數(shù)據(jù)。

[1] 謝菠蓀. 頭相關函數(shù)與虛擬聽覺[M]. 北京: 國防工業(yè)出版社, 2007.

XIE Bosun. Head related transfer function and virtual auditory [M]. Beijing: National Defence Industry Press, 2007.

[2] Fukudome K, Suetsugu T, Ueshin T. The fast measurement of head related impulse responses for all azimuthal directions using the continuous measurement method with a servo-swiveled chair[J]. Applied Acoustics, 2007, 68(8): 864-884.

[3] Kahana Y, Nelson P A, Petyt M. Boundary element simulation of HRTFs and sound fields produced by virtual sound imaging system[C]// AES 105 Convention, San Francisco, CA, USA, 1998.

[4] LIU Tiao, XIAO Qinghuo. Finite difference computation of head-related transfer function for human hearing[J]. J. Acoust. Soc. Am, 2003, 113(5): 2434-2441.

[5] Fels J, Buthmann P, Vorlander M. Head-related transfer functions of children[J]. Acta Acoustia United with Acustica, 2004, 90(5): 918-927.

[6] Wolfgang Kreuzer, Zhensheng Chen. A fast multipole boundary element method for calculating HRTFs[C]// AES 122ndConvention, 2007, May5-8.

[7] Lkulish Antani, Nico Galoppo, Adam Lake, et al. MIT database, The university of North Carolina at Chapel HILL, http://sound.media.mit.edu/resources/ KEMAR.html.

[8] “CIPIC database,” http://interface.cipic.ucdavis.edu.

[9] 仝欣, 齊娜. 橢球頭模與仿真頭模的指向性比較[J]. 電聲技術, 2012, 36(6): 43-46.

TONG Xin, QI Na. Comparison of the directivity between the ellipsoid head and the dummy head[J]. Audio Engineering, 2012, 36(1): 43-46.

[10] 祝家麟, 袁政強. 邊界元分析[M]. 北京: 科學出版社, 2009.

ZHU Jialin, YUAN Zhengqiang. Boundary element analysis[M]. Beijing: Science Press, 2009.

[11] Schenck H A. Improved integral formulation for acoustic radiation problems[J]. J. Acoust. Soc. Am, 1968, 44: 41-58.

[12] Przemyslaw Plaskota, Andrzej B. Dobrucki. Head-Related Transfer Function calculation using Boundary Element Method[C]// Audio Engineering Society, 2007, May 5-8.

[13] 杜功煥, 朱哲民, 龔秀芬, 等. 聲學基礎[M]. 南京: 南京大學出版社, 2001, 3.

DU Gonghuan, ZHU Zhemin, GONG Xiufen, et al. Acoustic foundation[M]. Nanjing: Nanjing University press, 2001, 3.

[14] 李增剛, 詹福良. Virtual.Lab Acoustics 聲學仿真計算高級應用實例[M]. 北京: 國防工 業(yè)出版社, 2010.

LI Zenggang, ZHAN Fuliang Virtual.Lab Acoustics acoustic simulation advanced application examples[M]. Beijing: National Defence Industry Press, 2010.

[15] “Structured Light Scanner”, http://www.david-laserscanner.com

[16] WU T W. Bonduary Elements in Acoustics[M]. WIT Press, Southampton, 2000.

[17] 333 Three D Systems Circle Rock Hill, SC 29730, USA, http://www.rapidform.com/

[18] Duda R O, Martens W L. Range dependence of the response of a spherical head model [J]. J. Acoust. Soc. Am., 1998, 104(5): 3048-3058.

Numerical method for calculating head-related transfer function of the Chinese dummy head BHead210

TANG Ling1, FU Zhong-hua2

(School of Computer Science, Northwestern Polytechnical University of China, Xi’an 710129,Shaanxi, China)

A numerical method for calculating the head-related transfer function (HRTF) of the Chinese dummy head BHead210is studied. With the Chinese dummy head BHead210, a structured light scanner is used to obtain its three dimensional data and then the boundary element method is used to calculate its HRTF. The calculation results with a rigid simple sphere model verify its correctness. The comparison results of the numerical calculation and those obtained by acoustic measurement are given. The Interaural Time Difference and Interaural Level Difference curves, different source-distance of HRTF amplitude spectrum are given too. The HRTF of the BHead21 is also compared to those of KEMAR reported by CIPIC.

head-related transfer function; Chinese dummy head BHead210; boundary element method

TN912

A

1000-3630(2014)-03-0237-06

10.3969/j.issn1000-3630.2014.03.011

2013-06-30;

2013-10-12

唐玲(1987－), 女, 海南儋州人, 黎族, 碩士研究生, 研究方向為 虛擬聽覺。

唐玲, E-mail: tangling923@mail.nwpu.edu.cn