郭 瑩， 閆美辰

(沈陽工業(yè)大學(xué) 信息科學(xué)與工程學(xué)院， 沈陽 110870)

基于鏡像源方法的室內(nèi)聲場脈沖響應(yīng)仿真*

郭 瑩， 閆美辰

(沈陽工業(yè)大學(xué) 信息科學(xué)與工程學(xué)院， 沈陽 110870)

針對基于室內(nèi)聲場脈沖響應(yīng)描述聲場音質(zhì)特性時其精確度不高的問題，以經(jīng)典鏡像源方法為核心，討論了基于變換頻域和反射系數(shù)的改進(jìn)算法.以矩形空間為具體實例，給定接收點和聲源位置，并利用計算機仿真室內(nèi)聲場脈沖響應(yīng).在不同參數(shù)條件下，對各個算法進(jìn)行了分析，估計了能量衰減曲線和房間混響時間.結(jié)果表明，基于變換頻域和反射系數(shù)的改進(jìn)算法具有更高的準(zhǔn)確性和可靠性.

室內(nèi)聲場； 脈沖響應(yīng)； 音質(zhì)特性； 鏡像源方法； 變換頻域； 反射系數(shù)； 能量衰減曲線； 混響時間

由于封閉聲場內(nèi)聲能密度、聲壓等參數(shù)都滿足齊次性和疊加性原理，因此，可將封閉的待測室內(nèi)聲學(xué)空間看作一個線性時不變系統(tǒng)[1].時域的通道特性可用聲源至接收點的聲場脈沖響應(yīng)來描述[2].借助現(xiàn)代信號處理技術(shù)，通過對脈沖響應(yīng)進(jìn)行測量和分析，可以得到聲學(xué)空間的各種聲學(xué)特性參數(shù)[3-5]，如混響時間、明晰度和聲場強度等，從而準(zhǔn)確地描述室內(nèi)聲場的變化規(guī)律.

當(dāng)仿真計算聲場脈沖響應(yīng)時，需要選取合適的聲場模型.鏡像源方法是以幾何聲學(xué)模型為理論基礎(chǔ)的經(jīng)典算法，可以較準(zhǔn)確地模擬聲場脈沖響應(yīng).自1979年以來，憑借理論簡單、易于理解和實現(xiàn)的優(yōu)點，鏡像源方法已被視為聲學(xué)研究的基礎(chǔ)理論，因而在聲學(xué)、工程學(xué)等很多領(lǐng)域中得到了廣泛應(yīng)用[6].可以利用鏡像源方法對封閉環(huán)境中聲音的傳播、封閉空間中噪聲的控制，以及聲場中的語言清晰度和語音傳輸指數(shù)進(jìn)行分析[7]，且鏡像源方法也可應(yīng)用于雙耳可視化和交互系統(tǒng)等設(shè)備的研發(fā)領(lǐng)域[8].

本文以經(jīng)典鏡像源方法為基礎(chǔ)，結(jié)合統(tǒng)計聲學(xué)，討論了幾種聲場脈沖響應(yīng)的具體算法，編制了仿真程序，給出了矩形空間聲場的計算實例，并對直接影響模擬結(jié)果的幾個重要方面進(jìn)行了討論.基于仿真得到了室內(nèi)聲場脈沖響應(yīng)，并估計了能量衰減曲線和室內(nèi)混響時間.

1 鏡像源方法

1.1 基本原理

圖1為鏡像源方法俯視圖.由圖1可見，在封閉空間中，由聲源S發(fā)射出來的聲波在空氣中傳播時，一部分聲波被接收點Q直接接收，另一部分聲波將被墻面或障礙物反射，經(jīng)過一定延時后再被Q接收.由于反射角等于入射角，故反射波的反向延長線與聲源到墻面的垂線之間必然存在某一交點Si(i=1，2，3，4)，反射聲波可以看做是從墻后的另一聲源Si發(fā)射出來的，因此，Si被稱為聲源S的鏡像源，且鏡像源Si與聲源S到墻壁的距離相等.當(dāng)獲知全部鏡像源的位置與能量(或聲壓)后，可以把聲源對于接收點的能量等效為相應(yīng)鏡像源的能量之和.

圖1 鏡像源方法俯視圖Fig.1 Vertical view of image source method

鏡像源數(shù)目與反射級數(shù)(聲源反射次數(shù))和壁面數(shù)目均有關(guān)系.隨著鏡像源數(shù)目的增加，總鏡像源數(shù)目將按照指數(shù)規(guī)律遞增.但是并非所有鏡像源都會對封閉聲場內(nèi)的接收點產(chǎn)生能量.事實上，產(chǎn)生能量的鏡像源只占很小的比例.因此，在鏡像源方法中，一個很重要的步驟就是判斷鏡像源的可見性，即要從全部的鏡像源中，確定哪一些鏡像源可以對接收位置產(chǎn)生能量.這一過程需要遵循的原則為產(chǎn)生鏡像源的壁面必須與鏡像源、接收點的連線存在實際交點.

圖2為鏡像源的可見性判斷示意圖，其中實線是聲音真正的傳播路徑，而虛線是所引輔助線.如圖2所示，對于聲源S和接收點Q而言，一級鏡像源S′是可見的，而二級鏡像源S″則是不可見的.通過鏡像源的可見性判斷，可以大大降低鏡像源方法的計算量.

圖2 鏡像源的可見性判斷Fig.2 Visibility judgment for image source

1.2 經(jīng)典算法

可以通過一系列可見鏡像源產(chǎn)生的聲壓和來模擬計算聲場脈沖響應(yīng).在矩形自由場中，接收點接收到的聲壓(頻域聲壓脈沖響應(yīng))可以表示為

(1)

根據(jù)鏡像源原理，鏡像源與接收點是關(guān)于某一墻面對稱的，因此，聲壓還可以表示為

exp(-jωt)

(2)

式中：R-為聲源到接收點的距離；R+為鏡像源到接收點的距離.此外，R-與R+可以分別表示為

(3)

(4)

考慮到房間中有6面墻，且每一個一階鏡像源還會產(chǎn)生高階鏡像源，計算比較復(fù)雜，可將式(2)改寫為

exp(-jωt)

(5)

式中：Rp為相對于坐標(biāo)原點而言，聲源到鏡像源的位置向量；Rr為考慮高階反射后的虛擬房間尺寸向量.此外，Rp和Rr可以分別表示為

Rp=(x±x′，y±y′，z±z′)

(6)

Rr=2(nLx，lLy，mLz)

(7)

Rp+Rr=(x±x′+2nLx，y±y′+2lLy，

z±z′+2mLz)

(8)

式中：Lx、Ly和Lz分別為房間的長、寬、高；n、l和m為相應(yīng)參數(shù).

對式(5)進(jìn)行傅里葉變換后，得到的時域脈沖響應(yīng)表達(dá)式為

(9)

式中，δ(·)為單位脈沖函數(shù).

在以上推導(dǎo)中，均需要假設(shè)墻面是堅硬的.在實際情況中，則需要考慮到每一個壁面均存在吸聲系數(shù)α和反射系數(shù)β的情況.根據(jù)Sabine公式[8]，可以得到二者的關(guān)系為

α=1-β2

(10)

當(dāng)考慮到房間反射系數(shù)后，可以將式(9)改寫為

(11)

式中：

Rp=(x-x′+2qx′，y-y′+2vy′，

z-z′+2kz′)

(12)

q、v和k為相應(yīng)參數(shù)，且分別等于0或1.

此外，室內(nèi)所有表面的聲學(xué)特性可用反射系數(shù)βx，ξ、 βy，ξ和βz，ξ來表征.其中：當(dāng)ξ=1時，表示緊鄰坐標(biāo)系原點的墻壁；ξ=2時，表示遠(yuǎn)離坐標(biāo)系原點的墻壁.

2 改進(jìn)算法

一般而言，由基本鏡像源方法仿真得到的脈沖響應(yīng)精確度較低，因而無法真實地反應(yīng)室內(nèi)脈沖響應(yīng).若把基本方法的室內(nèi)脈沖響應(yīng)表達(dá)式變換到頻域進(jìn)行計算，并將表達(dá)式中的反射系數(shù)定義為負(fù)反射系數(shù)，得到新的室內(nèi)脈沖響應(yīng)表達(dá)式，則會提高室內(nèi)脈沖響應(yīng)建模的精確度.

2.1 頻域計算

當(dāng)在時域處理離散信號時，很難把握時間點的選取.若選取離每個鏡像源傳播時間點最近的整數(shù)點，則時間點的選取較為粗略，因而會使脈沖響應(yīng)產(chǎn)生相當(dāng)大的誤差.當(dāng)在頻域處理信號時，信號本身是連續(xù)的，無需再把連續(xù)信號離散化，從而避免了時間點的選取問題.由式(11)可以得到聲壓脈沖響應(yīng)的頻域表達(dá)式，即

exp(-jωt)

(13)

此時，將式(13)進(jìn)行傅里葉逆變換后，可以得到時域狀態(tài)下室內(nèi)脈沖響應(yīng)表達(dá)式，即

hf(t)=F-1(Hf(ω))

(14)

式中，F(xiàn)-1為傅里葉逆變換.

2.2 負(fù)反射系數(shù)

設(shè)任何房間的表面吸收系數(shù)均由吸聲系數(shù)α表示，則由式(10)可以得到反射系數(shù)β的表達(dá)式[9-10]，即

(15)

由式(15)可知，β可以取正、負(fù)兩個值，在傳統(tǒng)鏡像源方法中，β值取的是正數(shù).然而，利用傳統(tǒng)鏡像源方法，得到的脈沖響應(yīng)與真實測量值相比，失真情況較為嚴(yán)重，無法清晰地反映室內(nèi)脈沖的真實情況.事實上，利用變換頻域與采用負(fù)反射系數(shù)的方法，仿真得到的脈沖響應(yīng)與實際情況更為接近.若將式(15)中的β定義為負(fù)數(shù)，可以實現(xiàn)對室內(nèi)脈沖響應(yīng)的精確建模.此時，可以得到新的聲壓脈沖響應(yīng)表達(dá)式，即

exp(-jωt)

(16)

此時，將式(16)進(jìn)行傅里葉逆變換后，得到的時域狀態(tài)下室內(nèi)脈沖響應(yīng)表達(dá)式為

hne(t)=F-1(Hne(ω))

(17)

3 矩形空間的算法驗證

矩形空間是最典型也是最簡單的一種三維模型，對其進(jìn)行研究和討論，有利于研究更為復(fù)雜的情形，因此，本文針對矩形房間進(jìn)行了仿真.假設(shè)矩形房間天花板材質(zhì)為石膏；地板材質(zhì)為樹脂.房間內(nèi)材料的吸聲系數(shù)如表1所示.

表1 房間內(nèi)材料的吸聲系數(shù)
Tab.1 Sound absorption coefficients of room materials

室內(nèi)空間吸聲系數(shù)四面?zhèn)葔? 90地面0 70頂棚0 99

以房間一角為原點，且以長為x軸、寬為y軸、高為z軸建立三維直角坐標(biāo)系，可以得到仿真室內(nèi)空間的相應(yīng)坐標(biāo).表2為仿真室內(nèi)空間的環(huán)境設(shè)置情況.

為了比較實驗結(jié)果的準(zhǔn)確性，在與實驗環(huán)境參數(shù)相同的房間中，獲取真實的室內(nèi)脈沖響應(yīng)，使之與測量結(jié)果進(jìn)行比較.首先采用幅頻特性以平坦譜的白噪聲信號作為激勵信號，然后記錄由揚聲器輸出的激勵信號在麥克風(fēng)接收端的脈沖響應(yīng)，最后通過解卷積的方法求得室內(nèi)脈沖響應(yīng).為了獲得較為準(zhǔn)確的脈沖響應(yīng)，可以通過多次記錄數(shù)值并取平均值的方式得到真實的室內(nèi)脈沖響應(yīng).

表2 仿真室內(nèi)空間的環(huán)境設(shè)置
Tab.2 Environment arrangement of simulated room space

空間尺寸聲源坐標(biāo)/m麥克風(fēng)坐標(biāo)/m采樣頻率/kHz聲速/(m·s-1)3 2m×4 0m×2 7m(2 0，3 0，2 0)(1 1，1 0，1 2)16343

圖3為由不同方法獲得的室內(nèi)脈沖響應(yīng)對比結(jié)果.由圖3可知，與真實的室內(nèi)脈沖響應(yīng)相比(見圖3c)，利用基本鏡像源方法獲得的室內(nèi)脈沖響應(yīng)的失真情況較為嚴(yán)重(見圖3a)；利用變換頻域與采用負(fù)反射系數(shù)方法得到的室內(nèi)脈沖響應(yīng)的精確度則相對較高(見圖3b).因此，利用變換頻域和采用負(fù)反射系數(shù)的鏡像源方法可以獲得更接近于真實室內(nèi)脈沖響應(yīng)的仿真結(jié)果.

為了進(jìn)一步比較經(jīng)過變換頻域和采用負(fù)反射系數(shù)的鏡像源方法與基本鏡像源方法的仿真準(zhǔn)確度，選取了5個不同位置的接收點進(jìn)行了分析，且各接收點的位置如表3所示.

表3 各接收點位置Tab.3 Locations of various receiving points m

為了進(jìn)行性能比較，采用歸一化誤差作為評價建模準(zhǔn)確性的準(zhǔn)則，且其表達(dá)式為

(18)

表4為由不同方法得到的室內(nèi)脈沖響應(yīng)誤差對比結(jié)果.

表4 室內(nèi)脈沖響應(yīng)誤差對比Tab.4 Comparison in errors of room impulse response dB

表4中改進(jìn)方法1是針對頻域方法的改進(jìn)；改進(jìn)方法2是針對頻域與反射系數(shù)的改進(jìn).由表4可見，無論接收點位置如何改變，與基本鏡像源方法相比，改進(jìn)鏡像源方法所得到的室內(nèi)脈沖響應(yīng)誤差更小，且改進(jìn)方法2獲得的室內(nèi)脈沖響應(yīng)更接近于真實情況.

4 應(yīng) 用

混響時間是廳、堂、館音響效果設(shè)計中一個重要的客觀評價參數(shù).利用脈沖響應(yīng)獲得聲場能量衰減曲線(EDC)后估計混響時間，已經(jīng)成為近年來國家標(biāo)準(zhǔn)化組織(ISO)對各聲學(xué)參數(shù)進(jìn)行標(biāo)定的標(biāo)準(zhǔn)測量方法.所謂能量衰減曲線是指聲音在充滿聲場后突然停止，某一點的聲能隨時間的推移逐漸衰減而形成的曲線.依據(jù)能量衰減曲線可以估計混響時間[9-10].

根據(jù)式(17)并結(jié)合Schroeder方法可以得到能量衰減曲線(EDC)的表達(dá)式

(19)

式中，E為能量.

根據(jù)式(19)可以仿真得到能量衰減曲線，結(jié)果如圖4所示.圖4中3條曲線是在不同的反射系數(shù)下得到的仿真結(jié)果，且反射系數(shù)的具體設(shè)置如表5所示.由圖4的能量衰減曲線可直接讀取出能量從-5 dB變化到-25 dB的時間.

圖4 能量衰減曲線Fig.4 Energy decay curve (EDC)

表5 室內(nèi)各表面的反射系數(shù)Tab.5 Reflection coefficients of various room surfaces

表6為分別采用了基于變換頻域和負(fù)反射系數(shù)的脈沖響應(yīng)方法、Eyring方法與Sabine方法[11]后獲得的混響時間數(shù)據(jù).表6中的誤差率是指估值與理想值作差取絕對值后，再與理想值作比后的數(shù)值.由表6可見，由脈沖響應(yīng)方法獲得的混響時間的準(zhǔn)確率遠(yuǎn)高于另外兩種方法.

表6 混響時間比較
Tab.6 Comparison in reverberation time

方法理想值/s估值/s誤差率/%脈沖響應(yīng)方法Eyring方法Sabine方法0 050 04840 100 082180 150 14440 050 045100 100 136360 150 212410 050 063260 100 138380 150 21543

5 結(jié) 論

本文以經(jīng)典鏡像源方法為基礎(chǔ)，討論了幾種聲場脈沖響應(yīng)的具體算法，編制了仿真程序，給出了矩形空間聲場的計算實例，并對直接影響模擬結(jié)果的幾個重要問題進(jìn)行了討論，最后基于仿真得到室內(nèi)聲場脈沖響應(yīng)，估計了能量衰減曲線和房間混響時間.仿真結(jié)果表明，利用基于變換頻域和負(fù)反射系數(shù)的鏡像源方法獲得的脈沖響應(yīng)的混響時間結(jié)果與理想值具有較好的一致性，且具有更高的準(zhǔn)確性和可靠性.

[1]成燁.室內(nèi)聲場的脈沖響應(yīng)分析 [D].太原：太原理工大學(xué)，2014.

(CHENG Ye.Impulse response analysis for room sound field [D].Taiyuan：Taiyuan University of Technology，2014.)

[2]馬宏偉，董明.基于矩形換能器空間脈沖響應(yīng)的相控陣次聲場研究 [J].機械工程學(xué)報，2014，50(18)：36-42.

(MA Hong-wei，DONG Ming.Research on acoustic field of ultrasonic linear phased array based on spatial impulse response of rectangular planar transducer [J].Journal of Mechanical Engineering，2014，50(18)：36-42.)

[3]黃坤朋，趙越喆.應(yīng)用時域有限差分法計算房間脈沖響應(yīng)和聲場參數(shù) [J].土木建筑與環(huán)境工程，2013，35(3)：157-161.

(HUANG Kun-peng，ZHAO Yue-zhe.Calculation of room impulse response and acoustic parameters by finite-difference time-domain method [J].Journal of Civil，Architectural and Environmental Engineering，2013，35(3)：157-161.)

[4]陳妍.長空間的聲學(xué)特性研究 [D].南京：南京大學(xué)，2011.

(CHEN Yan.A study on the acoustic of long enclosures [D].Nanjing：Nanjing University，2011.)

[5]曾向陽.室內(nèi)聲場計算機模擬發(fā)展40年(1968～2008) [J].電聲技術(shù)，2008(10)：15-17.

(ZENG Xiang-yang.Forty years of development of room acoustic computer simulation [J].Audio Engineering，2008(10)：15-17.)

[6]姜錦云.基于麥克風(fēng)陣列的近場和遠(yuǎn)場混合聲源定位 [D].成都：西南交通大學(xué)，2010.

(JIANG Jin-yun.Mixed near-filed and far-filed speech source local localization based on microphone arrays [D].Chengdu：Southwest Jiaotong University，2010.)

[7]焦偉奇.基于虛源法的聲矢量場建模與預(yù)報及軟件實現(xiàn) [D].哈爾濱：哈爾濱工程大學(xué)，2011.

(JIAO Wei-qi.Research on acoustic vector field calculation and analysis based on image source and software implementation [D].Harbin：Harbin Engineering University，2011.)

[8]矯姍姍.廳堂音質(zhì)控制計算機仿真的誤差分析 [D].青島：中國海洋大學(xué)，2013.

(JIAO Shan-shan.The error analysis of computer simu-lation for the hall acoustics quality control [D].Qingdao：Ocean University of China，2013.)

[9]Li K，Lam P.Prediction of reverberation time and speech transmission index in long enclosures [J].The Journal of the Acoustical Society of America，2005，117(6)：3716-3726.

[10]Omoto A.Measurement of sound absorption coefficients in a reverberation room [J].Journal of the Acoustical Society of Japan，2012：68(12)：616-621.

[11]劉延昌.廳堂聲場計算機模擬與實際測量結(jié)果對比分析研究 [D].邯鄲：河北工程大學(xué)，2011.

(LIU Yan-chang.Comparative analysis of the acoustic simulation and the filed measurement’s results [D].Handan：Hebei University of Engineering，2011.)

(責(zé)任編輯：尹淑英 英文審校：尹淑英)

Simulation for impulse response of room acoustic field based on image source method

GUO Ying, YAN Mei-chen

(School of Information Science and Engineering, Shenyang University of Technology, Shenyang 110870, China)

In order to solve the problem that the accuracy in describing the sound quality characteristics of acoustic field based on the impulse response of room acoustic field is not high, the classical image source method (ISM) was taken as the core, and the improved algorithms based on the transformation frequency domain and reflection coefficient were discussed. In addition, the rectangular space was taken as the specific example, the receiving points and sound source locations were given, and the impulse response of room acoustic field was simulated with the computer simulation. According to different parameters, the various algorithms were analyzed. In addition, the energy decay curve (EDC) and room reverberation time (RT) were estimated. The results show that the improved algorithm based on transform frequency domain and reflection coefficient has higher accuracy and reliability.

room acoustic field; impulse response; sound quality characteristic; image source method; transform frequency domain; reflection coefficient; energy decay curve; reverberation time

2015-12-17.

國家自然科學(xué)基金資助項目(61102123)； 遼寧省教育廳優(yōu)秀人才資助項目(LJQ2014011)； 遼寧省教育廳一般項目(L2014041)； 沈陽市科技局工業(yè)攻關(guān)項目(F14-073-2-00).

郭 瑩(1975-)，女，遼寧鐵嶺人，副教授，博士，主要從事非高斯信號處理與參數(shù)估計等方面的研究.

16∶06在中國知網(wǎng)優(yōu)先數(shù)字出版.

http：∥www.cnki.net/kcms/detail/21.1189.T.20160907.1606.012.html

10.7688/j.issn.1000-1646.2017.01.11

TN 912

A

1000-1646(2017)01-0055-06