相干聲場的快速預測方法研究*

呂巖，劉志紅，吳群，儀垂杰

（1.青島理工大學機械與汽車工程學院 青島，266520）

（2.工業流體節能與污染控制教育部重點實驗室 青島，266520）

引言

噪聲控制中，室內聲場的預測對于室內聲學響應及噪聲治理有重要意義。有效的聲場預估能獲得準確的聲場特征和主要噪聲源，為室內聲源的布置和制定合理的噪聲治理方案提供理論依據。

20世紀80年代，有學者提出眾多虛聲源模型來模擬室內聲場，這些預測模型忽略聲波之間的干涉效應，導致計算精度低、誤差大。為改善這一缺陷，Wang等［1］將分析不同聲波之間聲壓相干效應的虛源法應用于工廠車間和開放式辦公室的聲場預估，結果表明考慮聲波相干的模型能提高聲場預測精度，但因模型對聲源和邊界條件的簡化，誤差仍然突出。Lemire［2］建立了復雜虛源聲場預測模型，將聲波在邊界上的反射視為球面波。Nobile等［3］針對單次球面波反射，提出了數值求解方法。文獻［4-8］對不同聲學邊界條件的空間聲場進行預估，拓寬了復虛源模型的預測范圍，證明了將不同邊界簡化為不同的邊界阻抗并結合虛源法能有效預估聲場分布，減小預測誤差，但對大尺度空間、復雜聲源的預測，其計算精度有限且計算復雜度仍較高。

針對模型的預測精度，Wang等［9］對超大空間的聲場進行了研究，結果表明隨著空間體積的增大，聲場變化的主要原因是反射聲能量的改變，并基于虛源法與聲壓級衰減公式提出了超大空間聲場的預測公式，但超大空間邊界引起的散射更為復雜，使用不考慮聲波干涉的模型精度仍有限。徐祿文等［10］采用聲線法對影響變壓器主變室外聲場的因素進行了分析，其結果說明聲場邊界阻抗的分布對聲場的分布有重要的影響，準確的邊界條件是聲場精度的重要前提，但其提出的預測模型仍不考慮聲場的干涉。

為了提高預測精度，除了對預測模型進行改進，部分學者對聲源的處理也做了相應的研究。文獻［11-17］應用不同優化的噪聲源模型來重構分析設備的輻射聲場，預測結果的精度有了很大提升。等效噪聲源模型的優化研究從另一方面提高了聲場預估的精度，但同時復雜噪聲源模型在應用時的計算量也隨之增加。針對模型的計算效率，張朝金等［18］在射線模型和高斯射線束理論基礎上，利用并行化處理對Bellhop射線模型進行了優化并開發了BellhopMP模型，該模型通過多核計算機并行計算大大提高了水下聲場預估效率。此外，在聲學計算領域為提高計算效率，眾多學者將快速多極算法引入 計 算 之 中。Jiang等［19］將快速 多 極 算 法 引入2維聲學問題，并將計算復雜度由O(N2)降低到O(NlogN)。張炳榮等［20］將快速多極算法引入基本解法，實現了一種2維聲場快速預測方法。張士偉［21］將3維多極算法引入機械噪聲的預測，以發動機為例，說明該方法求解噪聲源輻射噪聲的有效性，但仍未見快速多極算法具體應用于相干聲場的預測當中。

針對封閉空間相干聲場預估保證精度同時滿足計算效率的要求，首先，基于ISM原理構建了虛擬接收點模型和考慮聲場相干性的基本預測模型，考慮聲場的多項因素以滿足精度要求；其次，引入快速多極展開算法，將多點對多點的復雜映射關系轉化為點集對點集的快速計算過程，對模型進行降階計算以提高計算效率；最后，通過實驗仿真驗證本研究方法在復雜聲場預測的高效性及適用性。

1 預測模型基本原理

1.1 復虛源相干模型

虛源法將噪聲源對室內接收點的作用轉化為若干階虛源對接收點作用總和，虛擬聲源與實際聲源關于邊界有一定的幾何對稱關系。虛源法原理如圖1所示，其中：S為聲源點；R為接收點；S(I)為虛擬源點。

圖1 虛源法原理圖Fig.1 Schematic diagram of the virtual source method

一個源點對接收點的作用相當于源點與無數虛源對接收點的作用之和，即

其中：P(x)為接收點聲壓；P0為直達聲貢獻；Pi為各虛源點貢獻。

直達聲的貢獻可直接由式P0=qG(rs，rr)計算，其中：q為源強；G(rs，rr)為格林函數。而虛源的作用考慮邊界的反射引入反射系數，即虛源的作用可表示為，Qj為邊界的單次反射系數，得到總的反射系數代入即可求得虛源對接收點的貢獻，依次求解所有虛源的貢獻代入式（1）可得接收點的聲壓。

筆者對于單次反射系數采用式（2）的復反射系數，即

其中：k為波數；βj為第j個墻壁的法向比聲導納；dN為該路徑到達接收點的傳播距離；θn為第j個界面上入射方向與界面法向的夾角；i為虛數單位。

I(k，dN，βj，θn)為關于k，dN，βj，θn的函數

根據不同邊界條件及傳遞路徑，由式（2）和式（3）計算出單次反射系數，進而求得單個虛源點與接收點的總反射系數，代入式（1）與直達聲貢獻相加即可求得接收點處的聲壓值為

式（4）為針對不同邊界條件的波疊加法基本原理，由此可知，噪聲設備在空間產生的聲場可由一系列配置的等效源疊加而成。等效源點的源強可由若干場點匹配而來，進而可用等效源點代替設備參與聲場計算。

根據式（1）計算各個等效源點對接收點的貢獻，求和即為所求聲場分布。但各源點與各自虛源點的可見性判斷、傳遞路徑分析及復反射系數的求解等已占有很大計算量，而等效源點的數目一般與等效精度有關，無法選取太少，因此對各等效源點采取虛源法的計算量是巨大的。鑒于此，引入快速多極基本算法，將各等效源點的作用簡化為一個參考點的貢獻，從而降低計算量。

1.2 快速多極展開算法

快速多極算法的基本思想是通過將核函數即本研究的格林函數多極展開、局部展開及相應的展開轉移，將源點集與接受點集中點與點的相互作用轉化為兩個點集單元的相互作用，原理如圖2所示。其中：yi代表噪聲點；xi代表接收點；X1，Y1為子晶格中心；X0，Y0為兩八叉樹結構的中心；M2M，L2L為子晶格中心間的作用關系；M2L為源點集與接收點集間的作用關系。現求解三維格林函數多極展開系數及展開系數轉移式，以下列出關鍵結論及推導式。

圖2 3維快速多極算法示意圖Fig.2 Schematic diagram of three-dimensional fast multipole algorithm

根據Graf加法原理，0階第1類漢克函數可展開為

其中：jn(·)為n階貝塞爾函數；()為n階1類漢克函數；Pn()為n次勒讓德多項式。

三維格林函數可表述為

其中：r″為xy間的距離；r′為Y1y間的距離；r為Y1x間的距離；γ為Y1x與Y1y的夾角。

根據上述，式（4）可表示為

另將格林函數在Y0處展開，得到展開系數轉移式

式（12）中，Wn，n′，m，m′，l為

通過多極展開系數式（9）及展開轉移系數式（10）即可將各子晶格作用轉移至父晶格，逐級由最底層轉移至最高層，從而實現將源點子集中各點對接收點的作用轉化為子集整體對接收點的作用。

2 相干聲場快速預測方法

本研究預測模型主要基于虛源法，將室內任意噪聲源轉化為等效源點集，引入快速多極算法計算虛源產生的聲場，而直達聲采用傳統基本解法。具體實現步驟如下。

1）對噪聲源設備進行等效處理。根據等效源思想，將設備產生的聲場等效為若干單極子點源的疊加聲場，若已知聲源邊界條件則用式（4）求解各源點源強；若未知則采用聲強或聲全息的方法重構聲源［16-17］。聲源的源強矩陣及位置矩陣分別為

其中：qi為第i個等效源源強；qi(x，y，z)為第i個等效源位置坐標。

2）對等效源點子集建立3維八叉樹結構。筆者主要針對源點子集進行簡化計算，只選取1個包含所有源強點的立方體，等分這個立方體得到第2層，依次對含有源點的小立方體進行劃分直至最末層立方體中的源點數量小于設定的閾值。根據八叉樹結構建立各級分層中心點的源強及位置坐標矩陣為

其中：yi，j為第i層第j個晶格中心的源強；yi，j(x，y，z)為中心的坐標。

根據等效源點所在晶格，將等效源點進行范圍劃分，一般晶格內源點數量閾值為1，即1個晶格內含1個等效源點。

3）建立虛源模型。若對等效源點子集建立虛源模型，則需建立相應階數的源點八叉樹結構，模型計算量較大，因此建立虛擬接收點模型，建立接收點關于邊界的虛擬接收點子集，而噪聲源對接收的作用表示為噪聲源對每個接收點貢獻的集合。在圖1中與虛擬源相似，可得虛擬接收點(xs，ys，zs)在空間的分布為

其中：(L，W，H)為邊界的尺寸。

為計算方便將坐標原點置于室內某一頂點，另l∈(-∞，∞)，n∈(-∞，∞)，m∈(-∞，∞)，因此給定了1組(l，n，m)的值就決定了1個特定的虛擬接收點。由此建立虛擬接收點的聲場響應值矩陣及位置矩陣為

其中：si為第i個虛擬接收點的響應值；si(x，y，z)為虛擬接收點位置。

由邊界的尺寸可得，邊界SL0，SL1，SW0，SW1，SH0，SH1在3維空間的虛擬邊界分布分別為SL0+2lL，SL1+2lL，SW0+2nW，SW1+2nW，SH0+2mH及SH1+mH，l，n，m∈(-∞，∞)。S分別表示各邊界面，當l，n，m=0時為6個實際邊界面。

4）上行歷遍。在源強八叉樹結構中，根據等效源點位置矩陣式（13）與八叉樹結構晶格中心位置矩陣式（15），利用式（9）與式（10）計算結構各層晶格的多極展開系數及轉移系數，將各等效源點作用逐級向父層轉移，直至最高層，將整個點集的作用轉化為中心Y0的作用QY0。

5）室內聲場計算。根據虛擬接收源位置矩陣與源強八叉樹結構中心Y0的坐標關系建立傳遞路徑函數，不同的虛擬邊界對應不同的邊界阻抗，求解路徑與空間虛擬邊界的交點從而確定該路徑的傳遞距離、反射次數、各單次碰撞反射系數與總的反射系數，建立各接收點對應的反射系數矩陣

逐次求解Y0對所有虛擬接收的貢獻

所有虛擬接收點均滿足多極算法的前提，根據步驟2和4建立的作用點Y0與步驟3建立的虛源模型計算所有反射聲的貢獻，根 據 步 驟1的等效源點模型用傳統基本解法計算直達聲P0，避免對多極算法的條件判斷。將直達聲與反射聲相加即得到室內任一場點的聲場信息。

3 實驗與仿真分析

3.1 仿真分析

3.1.1 仿真示例1

為對比FMA-ISM方法的計算精度，針對某長方體空間進行數值仿真模擬，仿真模型示意圖如圖3所示。為了便于對比其他方法，模擬長方形空間參數為1 000 mm×500 mm×500 mm，設置8個聲源和9個接收點，每個單極子聲源振幅A設為0.1，聲源及接收點的位置參數如表1和表2所示。墻面邊界阻抗屬性按照數據庫中普通廠房鋼架結構材料定義。

圖3 仿真模型示意圖Fig.3 Schematic diagram of the simulation model

表1 聲源位置參數表Tab.1 Location parameters of sound source mm

表2 接收點位置參數表Tab.2 Location parameters of receiver point mm

為驗證本研究方法的計算效率和精確度，現分別利用FMA-ISM方法、有限元方法及傳統聲線法對上述模型進行計算。由于基于波動理論的有限元方法精度較高，在此以有限元方法進行對比。FMA-ISM方法在Matlab下計算，計算取前5階虛源。聲線法利用Soundplan計算，Soundplan利用傳統的聲線法與虛源法，計算時不考慮聲波相干性。計算在同一設備進行，計算機CPUi5-8600，雙核6線程，運行內存為8 G。

本仿真算例中有限元網格劃分大小為10 mm，計算頻率范圍為100～5 000 Hz，計算步長為10 Hz，接收點R7計算結果對比如圖4所示。由圖可見，由于有限元方法與FMA-ISM方法均考慮相位影響，2種方法計算結果有相同的聲壓級變化趨勢，反映了聲波間的相位干涉。而有限元考慮了空間自身模態，其變化與FMA-ISM模型略有不同。另外，由于聲線法不考慮相位，只考慮能量的疊加，使得變化趨勢很小、誤差較大。

圖4 接收點R7計算結果對比圖Fig.4 Comparison of calculation results of receiving point R7

進一步對比3種方法，采用絕對誤差值對比FMA-ISM方法、傳統聲線法與有限元方法的誤差

其中：PFin為有限元計算結果；PMod為FMA-ISM方法計算結果；PSou為聲線法計算結果。

接收點R7誤差結果對比如圖5所示。由圖可見：FMA-ISM方法的絕對誤差基本在6 dB以下，個別與房間模態對應頻率達到7～9 dB，平均誤差為3 dB；聲線法的絕對誤差均在8 dB以上，部分達到了20 dB以上，平均誤差達到了12 dB。

圖5 接收點R7誤差結果對比圖Fig.5 Comparison of the error of the receiving point R7 result

表3為3種方法計算時間對比。由表可知，FMA-ISM方法在保持一定誤差精度的前提下，計算效率較有限元方法提高50%以上。但是本算例模型較為簡單，隨著模型復雜度的進一步提高，有限元方法的計算時間會大大增加。另外，聲線法計算時間雖然也較低，但結果誤差較大，不滿足計算精度要求。

表3 不同方法計算時間對比Tab.3 Comparison table of calculation time by different methods s

考慮在同一頻率（500和4 000 Hz）下的各接收點結果，如圖6和圖7所示。由圖可見，本研究方法計算結果均與有限元結果變化趨勢趨于一致，而聲線法計算結果變化趨勢不明顯且部分接收點差值較大。結果對比反映出聲波的相干效應對聲場的影響較大，尤其當聲源環境與空間環境均較復雜時，忽略相位影響將導致結果有誤。

圖6 接收點R1～R9在500 Hz計算結果對比圖Fig.6 Comparison of calculation results of receiving points R1～R9 at 500 Hz

圖7 接收點R1～R9在4 000 Hz計算結果對比圖Fig.7 Comparison of calculation results of receiving points R1～R9 at 4 000 Hz

考慮各方法的誤差因素。有限元方法的計算效率直接與空間大小和網格大小有關，同時仿真計算需要對邊界條件進行合理的簡化，且有限元聲源輸入條件難以準確獲取造成誤差。Soundplan作為一種大型聲場仿真軟件，固化的程序使得操作較為簡單，計算廠區大范圍的聲場計算效率高，但其精度受算法限制，達不到誤差要求。FMA-ISM方法建立在現有虛聲源的基礎上，計算代碼簡單，易實現，誤差因素同有限元相同，在準確的聲源與邊界條件下，能獲得較好的計算結果。

3.1.2 仿真示例2

為進一步對比FMA-ISM方法與傳統ISM方法的誤差及計算效率，另針對某矩形空間聲場進行仿真計算。仿真模型如圖8所示。房間尺寸參數為10 m×5 m×5 m，聲源隨機分布在x∈(2 m，3 m)，y∈(2 m，3 m)，z∈(2 m，3 m)的立方體內，聲源由單極子源構成，單個源強100 dB，數目分別為8，16，32，64，128及256，接收點位于(7.5 m，2.5 m，2.5 m)。為對比空間參數對模型計算精度影響，現分別取3類空間邊界：①空間各邊界均為理想剛性邊界；②空間各邊界吸聲系數均為0.1；③空間頂部吸聲系數為0.3，其余邊界均為0.1，各邊界的吸聲系數處處相等。計算在同一設備進行，計算2 000 Hz聲場分布，計算過程取5階虛源，計算結果與誤差對比如圖9和圖10所示。其中：黑色線為ISM方法計算結果；紅色線為FMA-ISM方法計算結果。

圖8 仿真模型示意圖Fig.8 Schematic diagram of the simulation model

圖9 仿真計算結果對比Fig.9 Comparison of simulation results

圖10 仿真計算結果誤差對比Fig.10 Error comparison of simulation results

根據計算結果，在剛性邊界下FMA-ISM方法與傳統ISM方法結果趨于一致，絕對誤差在2 dB內。隨著邊界條件的復雜，兩者計算誤差增大，第3類邊界條件計算結果最大誤差達6 dB。分析原因為FMA-ISM方法簡化源點的作用為源點集合的作用，針對聲源的可見性判斷、總吸聲系數的求解等僅在源點集中心計算，而傳統ISM方法針對每一源點逐個分析計算，由此產生誤差，且當邊界條件差別越大，誤差越大，但二者平均誤差基本滿足工程計算要求。

各類邊界的計算時間如表4所示。由表可知，傳統ISM方法因對源點逐個分析，導致計算時間隨著源點數目的增加急劇增加；FMA-ISM方法隨源點數目的增加，計算時間增量較少，增加時間主要為源點集八叉樹結構復雜度的增加，而傳播分析時間增加較少。因此，改進的FMA-ISM方法可有效提高ISM方法的計算效率。

表4 不同模型計算時間對比表Tab.4 Comparison table of calculation time of different models s

3.2 實驗分析

為驗證FMA-ISM方法對聲場預估的適用性，在消聲室針對某噪聲設備進行聲場預估分析。消聲室邊界條件簡單，易實現計算要求，實驗測試現場如圖11所示，所用消聲室為金屬尖劈半消聲室，消聲室尺寸為10 m×8 m×3.7 m。噪聲源設備為某工業除濕機，除濕機尺寸為0.4 m×0.5 m×1.5 m。采用聲強法對除濕機本體噪聲進行測試，除濕機聲功率為94.77 dB。

圖11 實驗測試現場圖Fig.11 Experimental test scene

將除濕機置于消聲室頂角處，測量除濕機水平5 m×5 m范圍內聲場分布及東側2 m處2 m×2 m垂直立面的聲場分布。水平測點高度為1.2 m，測點間距為0.3～0.4 m，共200個測點，測點布置如圖12所示。垂直立面測點間距為0.2 m，共100個測點。測試用儀器為Norsonic150聲振分析儀，聲強探頭為GRAS 50AI-C。水平面和垂直立面的聲場分布分別如圖13和圖14所示。

圖12 現場測點布置圖Fig.12 Layout of on-site measuring points

圖13 水平面聲場分布圖Fig.13 Horizontal sound field distribution diagram

圖14 垂直立面聲場分布圖Fig.14 Distribution diagram of vertical elevation sound field

根據測試結果，在消聲室內四周及頂部邊界的吸聲系數較高，反射聲主要由底部反射貢獻，聲場分布與自由場接近。運用等效源基本思想，將除濕機等效為規則布置的單極子點源組合，等效源點數量分別為8和64。分別利用FMA-ISM及ISM方法對立面的聲場進行預估計算，選取場點中35個測點作為評價點，實測值與預估結果如圖15所示。圖中A類折線為現場實測值；B，D類折線為FMA-ISM方法，等效源點數分別為64和8的預估計算值；C，E類折線為ISM方法，等效源點數分別為64和8的預估計算值。根據對比結果，隨著等效源點數的增加，預測值與實測值趨于一致，ISM方法的計算結果精度高于FMA-ISM方法。為進一步分析誤差，采用絕對誤差對比分析各模型預估精度

圖15 選點計算結果對比圖Fig.15 Comparison of calculation results of point selection

其中：Pcal為理論計算值；Pre為實際測量值。

圖16為誤差率對比圖。圖中A，C類折線為FMA-ISM方法，源點數為64和8的預測誤差；B，D類折線為ISM方法，源點數為64和8的預測誤差。當源點數增加到64時，FMA-ISM方法與ISM方法的誤差均保持在2 dB內。隨著等效源點數量的增加，模型誤差減小，源點數增加到一定量后模型誤差趨于一定值，但隨著源點數的增加，模型計算量也隨之增加。表5為2類模型的計算時間，可以看出，隨著源點數的增加，ISM方法的計算時間急劇增加，而FMA-ISM方法的增加量較小，這與仿真示例結果一致。本研究實驗所用除濕機尺寸較小，源點數增加至64已可滿足要求。針對不同噪聲源設備，根據設備尺寸及噪聲產生機理合理選取等效源點數目，既可保證模型計算的精度也使得計算量控制在一定范圍內。通過實驗證明FMA-ISM模型可應用于現實的聲場預估。

圖16 選點計算結果誤差對比圖Fig.16 Comparison of errors in calculation results of point selection

表5 不同模型計算時間對比表Tab.5 Comparison table of calculation time of different models s

4 結論

1）通過復虛源方法與快速多極算法構建了封閉空間的聲場預估模型（FMA-ISM），以解決封閉空間聲場預估精確度不高、計算效率低等問題。通過復虛源原理與等效源方法建立虛擬接收點模型，對聲源及邊界條件進行等效后預估聲場，預估結果與實驗實測值誤差保持在3 dB內，能滿足預估精度的要求。引入快速多極算法對模型進行簡化計算，與其他方法對比，針對一定復雜度的空間聲場計算效率可提高50%以上，在保持預估精度的前提下大大提高了計算效率。

2）采用實驗與仿真對模型進行驗證，并與其他預估模型及方法進行對比。仿真1驗證了FMA-ISM方法結果與精度較高的有限元法結果趨于一致，兩者平均誤差在3 dB內，精度可滿足工程應用要求。仿真2驗證了在相同計算條件下，FMA-ISM通過對虛源貢獻的降階分析，計算效率較傳統ISM方法提高50%以上。最后，通過實驗室示例驗證了FMA-ISM方法在室內聲場預測的高效性及適用性。