一種基于麥克風陣列用于分離單極子和偶極子聲源的方法

周緯,楊明綏,馬威

1.上海交通大學 機械與動力工程學院,上海 200240

2.中國航空發動機集團有限公司 沈陽發動機研究所,沈陽 110015

3.上海交通大學 航空航天學院,上海 200240

4.上海交通大學 燃氣輪機與民用航空發動機教育部工程研究中心,上海 200240

隨著社會的快速發展,噪聲對人類生活、工業生產造成的弊端日益受到重視,如何更好地識別潛在的噪聲源位置及其強度,對于降噪設計尤為重要?；谘舆t疊加的波束成形算法CB(Conventional Beamforming),具有簡便、計算速度快和魯棒性高的優點,被廣泛應用于噪聲源的位置和強度估計。但是在實際工程應用中,噪聲源可能同時包含單極子源和偶極子源,應用單極子傳播假設的波束成形算法求解時,無法準確進行聲源定位。目前學術界還沒有有效的聲源定位技術,可以分離單極子源和偶極子源。

為了有效地對組合聲源進行分離,可以從反卷積算法和聲源傳播模型入手。由于受到陣列孔徑的限制,波束成形算法在低頻情況下的空間分辨率較低。因此,具有高分辨率的反卷積算法,如DAMAS (Deconvolution Approach for the Mapping of Acoustic Sources)算法、CLEAN算法,及其拓展算法DAMAS-C算法、CLEANSC(CLEAN for Spatial Source Coherence)算法應運而生。在航空聲學應用中,美國國家航空航天局(NASA)的Langley噴流噪聲實驗室成功運用DAMAS算法對幾種不同構型的噴嘴噴流噪聲進行了測量,比其他反卷積算法取得了更好的效果,此后DAMAS算法被廣泛地運用于航空聲學研究。由于DAMAS 具有迭代計算量大、計算成本高的缺點,Ma和Liu通過壓縮網格提高了DAMAS算法的計算效率。

另一方面,波束成形算法是通過建立聲源傳播模型,根據聲源輻射在空間傳播函數推導出聲源位置,所以要求定位聲源模型與所建立的求解模型要對應,而構建聲源模型往往受聲源類型、運動狀態影響。通常導向矢量采用的是基于單極子傳播特性的假設,因此采用標準波束成形技術無法準確估計一個簡單偶極子的位置和強度,Jordan等發現了這個問題,并對信號處理過程進行修正以考慮偶極子的傳播特性。為了修正偶極子假設模型,Liu等進一步推導了偶極子的聲源輻射傳播模型,提出了一種新的偶極子波束成形,改進了以前的信號校正方法,并利用實驗和仿真驗證了該方法恢復偶極子源的能力。此外,如果聲源處在運動中,由于存在多普勒效應,用于靜止聲源定位的反卷積算法無法有效識別運動聲源,因此產生了運動波束成形算法。

在實際應用中,聲源可能由不同類型、不同運動狀態的子聲源組成,因此學術界尋求可以分離組合聲源的成像算法。為了處理由靜止聲源和運動聲源組成的組合聲源,Mo 和Jiang基于DAMAS算法,在波束成形輸出和實際聲源分布之間構造一個線性方程組,成功將一個靜止源與一個運動聲源分離開來。當聲源數目較多時,Schmidt和Franks在傳統波束成形的基礎上,將導向矢量用導向矩陣表示,矩陣的每列對應不同模型假設的聲源。在涉及聲源相干性的問題時,Yardibi等提出一種用于相干聲源定位的協方差擬合方法,可以有效計算復雜度較低的非相干、部分相干乃至相干聲源。Suzuki進一步設計了一種廣義逆波束成形算法,通過將交叉互譜矩陣分解成特征模提取每個相干信號,利用廣義逆技術給出每個聲源的參考解。受Suzuki啟發,Pan等也是利用交叉譜矩陣分解,設計了一種多極子正交波成形算法,通過對特征向量的歸一化消除了聲源的指向性,用球面諧波表示聲源,聲源的類型和方向可以從重建的輻射模式中推導出來,可以更準確地估計聲源的位置、類型和方向。

盡管學術界在識別多種聲源上做出了很多努力,但是還沒有比較有效的成像算法可以將組合聲源分離開。受Mo和Jiang分離運動和靜止聲源的啟發,本文提出一種混合反卷積算法Hybrid DAMAS(HD),用于分離包含單極子源和偶極子源的組合聲源。所設計的混合反卷積算法建立在DAMAS算法上,并采用壓縮網格法來優化計算效率。在實際聲源分布與波束成形輸出間構建一個線性方程組,通過求解線性方程組,同時從單極子波束成形和偶極子波束成形中提取單極子源和偶極子源。本文設計了4個數值仿真算例及2組實驗來檢驗方法的有效性,考察定位位置誤差和強度誤差2個指標,并與波束成形算法和DAMAS算法作比較。

1 理論分析

1.1 聲源模型

給定一個聲源在具有穩定流動的均勻介質中位置ξ(1,2,…,),直角坐標系下的麥克風位置x (1,2,…,),馬赫數列向量/,其中為流體流速。計算從聲源到采集麥克風的傳播矢量R=x -ξ,設置1、聲源頻率、聲速。

計算單極子聲源在流體中的傳播時間Δ,以及在直角坐標系下的格林函數g

式中：為聲源角頻率。

一個偶極子聲源可以認為是由2個相位相反的單極子聲源組成,設置2個單極子源距離矢量,對于給定的,求解頻域內偶極子近場波動方程

式中：()為偶極子的聲壓幅值；Δ(·)為拉普拉斯算子。

對式(3)進行簡化,得

則在直角坐標系下,有穩定流動的偶極子格林函數表示為

1.2 混合反卷積算法

提出的混合反卷積方法需要在波束成形輸出和實際聲源分布之間構建一個線性方程組,構建4個點擴散函數矩陣PSF (Point Spread Function)：為偶極子波束成形和偶極子源之間的PSF矩陣,為偶極子波束成形和單極子源之間的PSF矩陣,為單極子波束成形和偶極子源之間的PSF矩陣,為單極子波束成形和單極子源之間的PSF矩陣。

計算傳播矢量矩陣：

將傳播矢量矩陣中的列向量進行歸一化,得到導向矢量的列向量

則導向矢量矩陣表示為

掃描平面上的聲壓分布用表示,假設聲源互相獨立,在掃描網格點中設置偶極子聲源與單極子聲源,其余網格點聲壓設置為0,為

根據傳播矩陣,聲源信息被傳遞給了麥克風采集系統,并通過計算機記錄成時域信號,接著利用FFT(Fast Fourier Transform)將時域信號轉換成頻域信號,最后利用式(14)計算出交叉互譜矩陣

傳統波束成形輸出結果：

式中：上標H 為轉置共軛。

將式(13)、式(14)代入式(15),波束成形輸出結果可以改寫為

式中：為源強度。傳播矢量矩陣攜帶聲學模型信息,導向矢量矩陣攜帶波束成形信息。

計算混合點擴散矩陣

其中：

式中：為偶極子傳播矢量矩陣；為單極子傳播矢量矩陣。

構造模擬聲源聲壓分布與波束成形輸出結果之間的線性方程組：

將式(18)代入式(19),混合反卷積線性方程組可以表示為

式中：［］為偶極子源和單極子源聲壓分布的列向量；［］為偶極子波束成形和單極子波束成形輸出結果的列向量。

2 仿真設置

設計4個仿真算例,用于驗證本文提出的混合反卷積分離算法的有效性。如圖1 所示,算例1～算例3設置了由一個單極子源和一個偶極子源組成的雙點源,2個聲源強度差異分別設置為0、10、-10 d B。算例4考察在多聲源狀況下該方法的有效性,設置了8個聲源,聲源強度均為94 d B。假設點源位于掃描平面上,如圖1所示,假設單極子源和偶極子源對于組合聲源的影響是獨立的。掃描平面大小為0.6 m×0.6 m,網格尺寸為0.01 m×0.01 m,共計有3 721個網格點。采用具有40個通道的螺旋麥克風陣列記錄壓力信號。測量平面與聲源掃描平面平行,并放置在距離掃描平面0.6 m 的位置,陣列中心位于軸上。采用DAMAS算法迭代1 000次,并采用文獻[14]中的壓縮網格方法來提高計算效率。聲源分布圖所展示的動態范圍為20 d B。圖1顯示了在均勻介質中流動的(0.5,0,0)和偶極子源距離矢量(0,0.02,0)m,所有的算例都在5 k Hz頻率下進行研究。

圖1 數值仿真聲源分布Fig.1 Numerical simulation sources distribution

所使用的仿真點源的聲源類型、位置和強度見表1、表2。

表1 雙點源分布信息Table 1 Double point sources distribution

表2 算例4多點源分布仿真信息Table 2 Multi-point sources distribution for Case 4

2.1 雙點源等強度分布

算例1中,偶極子源S1與單極子源S2等強度為94 dB,圖2展示了不同算法下算例1的聲源分布云圖,云圖中反映強度分布的物理量為聲壓級(,參考聲壓取2×10Pa)。圖2(a)表示基于單極子假設的波束成形求解結果,單極子源被準確識別出來,偶極子源被抑制,并沒有出現。在圖2(b)中,可以比較清晰地觀察到偶極子的存在,聲源強度也基本準確,但是仍然可以在云圖上看到單極子源。波束成形的空間分辨率如往常一樣比較糟糕,存在大量的偽源。圖2(c)、圖2(d)利用基于偶極子源和單極子源假設的DAMAS 算法分離組合源,提高云圖的分辨率。即使在DAMAS提高了結果云圖的空間分辨率后,仍然在圖上存在較多的旁瓣。偶極子源和單極子源分別作為干擾噪聲出現在對方的云圖上,聲源的強度誤差較小。圖2(e)、圖2(f)采用本文HD 算法分離組合源。圖2(e)展示的單極子聲源分布結果,比較準確地分離出了單極子源。圖2(f)展示了通過本文HD 算法分離出的偶極子聲源分布結果,在聲源強度上比較接近預設值94 d B,但在聲源位置上存在一定偏差,位置誤差為0.01 m。算例1仿真結果表明HD 算法成功地分離了2種等強度分布的不同類型的聲源。

圖2 算例1聲源定位云圖(f=5 k Hz)Fig.2 Sound source localization contour map for Case 1(f=5 k Hz)

2.2 偶極子聲源強度占優,S1強度比S2大10 dB

在算例2中,偶極子源的強度比單極子源大10 dB,考察當偶極子源比單極子源強時本文方法的有效性。圖3(a)偶極子源基本被過濾掉了,單極子比較清晰地出現在預設的聲源位置處,但是周圍有比較明顯的旁瓣,強度符合預設的84 d B。在由偶極子波束成形計算出的結果云圖中,如圖3(b)所示最大聲源出現在偶極子處,并且同時對單極子源有了很好的抑制作用。DAMAS算法很好地提高了云圖的分辨率,在圖3(c)中,單極子DAMAS消除了偶極子的旁瓣干擾保持源強,而單極子卻被當作偶極子處理并且伴有旁瓣；在圖3(d)中,盡管偶極子源作為最大聲源被準確識別,但少量的偽源仍然存在分布云圖上。本文HD 算法的計算結果如圖3(e)、圖3(f)所示,在偶極子聲源強大比單極子聲源強度大的情況下,它可以成功地區分2種不同類型的源,同時保持聲源強度的準確性。在位置誤差上,圖3(f)中偶極子位置存在0.01 m 的偏差。

圖3 算例2聲源定位云圖(f=5 k Hz)Fig.3 Sound source localization contour map for Case 2(f=5 k Hz)

2.3 單極子聲源強度占優,S1強度比S2小10 dB

算例3與算例2情況相反,偶極子的聲源強度比單極子小10 d B,考察當偶極子源比單極子源弱時該方法的有效性。圖4(a)中,單極子較完美地從組合聲源中被提取出來,保證了定位位置和強度的準確性。圖4(b)中,偶極子源被識別為最大聲源,周圍伴隨著較多的旁瓣,強度比預設值大6 d B,單極子附近也存在比較多的旁瓣。從圖4(c)、圖4(d)可以發現,基于偶極子假設的DAMAS無法有效分離出偶極子聲源,而基于單極子假設的DAMAS可以成功分離出單極子聲源。利用本文HD 算法分離該算例下的組合聲源,成功地將單極子源和偶極子源從組合源中分離開,保持定位強度的準確性展示在各自的潔凈聲源分布圖上,偶極子源的位置仍然有0.01 m的偏差,如圖4(e)、圖4(f)所示。

圖4 算例3聲源定位云圖(f=5 k Hz)Fig.4 Sound source localization contour map for Case 3(f=5 k Hz)

2.4 多個聲源

算例4設置8個聲源,包含4個偶極子源和4個單極子源,考察本文算法在處理多源問題時的有效性。組合聲源的位置分布圖和類型如圖1(b)、表2所示,8個聲源強度均為94 dB。在圖5(a)中,4個等強度的單極子源可以清晰地被識別出,盡管分辨率較差,但還是很好地過濾了偶極子聲源。圖5(b)展示了4個偶極子源的位置,整個云圖分辨率較差,存在比較多的旁瓣和偽源,有些位置單極子源作為干擾存在。DAMAS算法求解結果如圖5(c)、圖5(d)所示,盡管DAMAS提高了聲源分布圖的空間分辨率,但依然存在大量的偽源,所有的聲源定位都比較混亂,整個云圖結果都比較差。運用本文HD 算法對這8個聲源進行分離,預設的點源都可以在圖中找到,圖5(e)顯示了單極子源的分布,可以清楚地觀察到聲源的強度相差不大,但是在位置上與預設位置存在一定的偏離。圖5(f)中的偶極子源分布也比較完美,可以準確清晰地觀察預設的4個偶極子,定位位置上還是存在一些偏差。

圖5 算例4聲源定位云圖(f=5 k Hz)Fig.5 Sound source localization contour map for Case 4(f=5 k Hz)

3 實驗設置

為了進一步驗證本文算法,在上海交通大學航空發動機研究院的消音風洞上進行實驗驗證。通過在噴嘴出口外設置一根平行出口截面放置的均勻圓柱,利用圓柱擾流生成偶極子,并通過放置一個音響作為單極子,從而構造出包含單極子和偶極子的組合聲源。實驗中的噴嘴出口尺寸25 cm×25 cm,使用的空心不銹鋼圓柱直徑為6 mm,放置在距離出口截面20 cm 處。采用一個具有64通道的螺旋麥克風陣列,陣列平面與噴流中心線平行,距離1 m,使用了38個麥克風傳感器,型號為BK-4958,布置在螺旋陣列前38個通道上。數據采集機箱型號為NI PXle-1082,機箱主板型號為NI PXle-8840,采集板卡為NI PXle-4499,采樣頻率48 k Hz,持續10 s,頻譜計算方法采用Welch法,FFT 分塊長度為1 024。掃描平面尺寸1 m×1 m,網格大小為0.02 m×0.02 m,聲源分布圖動態范圍為10 d B。

3.1 實驗布置及頻譜分析

實驗布置如圖6(a)、圖6(b)所示,分別將不銹鋼圓柱橫向、縱向布置,依次產生一個向偶極子源和一個向偶極子源,在圓柱中心上方40 cm 處布置了一個音響作為單極子源。首先使用NI(National Instruments)采集器采集環境噪聲,接著打開藍牙音箱播放風扇葉片白噪聲,采集單極子的信號,然后啟動消音風洞,調整來流速度在110～115 m/s,采集組合聲源的信號,關閉音箱,采集偶極子信號,最后逐步減小流速,關停風洞。

圖6 現場聲源分布Fig.6 Sound source distribution

分別對以上2 種布置進行頻譜分析,由于FFT 的結果是復數,伴隨著相位,因此對于38個麥克風通道不能進行簡單的線性平均,這里采用的是幾何平均,對個通道的積求次方根,幾何平均頻譜如圖7、圖8所示。

圖7 實驗1頻譜圖Fig.7 Spectrum for Experiment 1

圖8 實驗2頻譜圖Fig.8 Spectrum for Experiment 2

不難發現,偶極子聲源在組合聲源中占優,已知圓柱直徑、氣流速度,根據斯特勞哈爾數可以推算出圓柱繞流產生的偶極子發聲頻率在3 200 Hz左右,與頻譜圖大致吻合,此外在實際測試過程中還發現在3、6、9 k Hz下存在較明顯的尖峰純音凸起,應是風洞變頻器干擾噪聲,在1 k Hz附近也出現了一個較大的干擾,應是本實驗中采用的風洞噴嘴出口以及夾具與氣流產生的干擾噪聲。下面對采樣信號進行聲源定位,分析頻率選取采樣信號聲源強度最大的頻段1 k Hz,考慮到波束成形算法在低頻下的低分辨率與混疊現象,又選取寬頻5 k Hz下的結果進行對比。

3.2 聲源定位云圖

圖9～圖12 分別為選取中心頻率1、5 k Hz時對2種偶極子布置的聲源定位結果。通過對頻譜分析可知,信號的最大聲源在1 k Hz左右,因此選取了1 k Hz作為分析頻率。首先看圓柱橫向布置,即偶極子的指向性為軸,CB 對偶極子進行了分離,盡管存在分辨率較差的情況,而對于單極子的識別中,未能消除偶極子對單極子定位的干擾。DAMAS 對于組合聲源進行了有效分離,在各自的定位云圖上都消除了對方的干擾,但是定位出來的聲源位置有較大的誤差,并且存在偽源。HD算法對DAMAS的結果進行了優化,聲源的定位位置精度有了提高,并且旁瓣、偽源抑制能力有了改善。對于圓柱縱向布置的情況,即偶極子指向性為軸,由于圓柱豎放時,氣流打在桿上形成了一條線性聲源,聲源強度可能超過單極子聲源強度,在此分析頻率CB、DAMAS以及HD都沒有能夠有效分離組合聲源。

圖9 實驗1聲源定位云圖(f=1 k Hz)Fig.9 Sound source localization map for Experiment 1(f=1 k Hz)

圖12 實驗2聲源定位云圖(f=5 k Hz)Fig.12 Sound source localization map for Experiment 2(f=5 k Hz)

對于寬頻5 k Hz的分析頻率,當圓柱橫向放置時,隨著頻率的升高,各個算法的分辨率有了顯著的提升。CB比較清晰地定位了偶極子,但是并未消除單極子的干擾,而在單極子的定位中,由于偶極子的強度較大,將偶極子誤判為最大的聲源。DAMAS 提升了CB 的分辨率和動態范圍,但是無論在識別單極子還是識別偶極子都未能削弱對方的干擾,2個聲源都比較清晰地被定位出。HD 改善了DAMAS的結果,在定位偶極子的同時盡可能的消除了單極子,而對于識別單極子中存在的較強偶極子也進行了一定的削弱。當圓柱豎直放置時,由于聲源與陣列的相對位置影響,CB將偶極子識別成了一個線性聲源,對于偶極子和單極子的定位相差不大。DAMAS反映的是CB結果在分辨率與動態范圍的改善,但是對于次聲源的消除表現一般,在偶極子與單極子云圖上都可以清晰地看見2個聲源。HD 算法進一步優化了DAMAS結果,準確定位偶極子的同時消除了單極子,盡管在定位單極子時沒有能夠完全消除偶極子,但是相較于DAMAS對單極子的定位結果已經有了改善。

圖10 實驗2聲源定位云圖(f=1 k Hz)Fig.10 Sound source localization map for Experiment 2(f=1 k Hz)

圖11 實驗1聲源定位云圖(f=5 k Hz)Fig.11 Sound source localization map for Experiment 1(f=5 k Hz)

4 討 論

基于DAMAS設計了一種混合反卷積算法,通過螺旋麥克風陣列測量將單極子源和偶極子源相結合的聲源分離開來。在對聲源類型的假設處理中,還存在以下問題：

1)偶極子指向性。本文方法側重于分離單極子聲源和已知指向性的偶極子聲源,例如風扇旋轉造成的旋轉偶極子聲源,如何分離單極子和未知指向性的偶極子聲源,將是未來的一個研究方向。

2)聲源的相干性。本文方法是基于假設單極子聲源和偶極子聲源相互獨立的,尚不能處理存在較強相互干涉的組合聲源。未來本文方法有望通過建立四極子模型應用于噴流噪聲的成分識別,但是由于真實噴流中各聲源之間存在較強大的相互干渉,因此將來有必要對如何分離具有相互干涉的聲源進行深入研究。

5 結 論

提出了可以實現單極子源和偶極子源分離的方法,這是當前基于單一聲源類型假設的波束成形算法和DAMAS 算法所不能做到的。本文方法建立在各類型聲源相互獨立的基礎上,在構造線性方程組的時候,同時考慮了單極子源和偶極子源的輸出,通過求解該方程組的解得到了單極子和偶極子的聲源分布。仿真和實驗證明,本文方法不僅可以在原始網格上顯示不同類型聲源的分布,而且可以保持較低的定位位置偏差和強度估計誤差。即使在處理較多聲源時,也能較好地區分開不同類型的聲源并且保持較明顯的源強差異。盡管混合反卷積算法被證明是分離由單極子和偶極子組成的組合源的有效方法,但在定位偶極子源過程中仍然存在很小的偏差。

本文方法有望應用于航空發動機氣動噪聲識別。航空發動機的噴流噪聲常常由單極子、偶極子及四極子源組成,通過在混合反卷積中添加四極子聲源模型,有望從噴流噪聲提取目標源,以加深對噴流噪聲產生機理以及后續降噪設計的研究。