似P范數特征值分解高分辨率聲源定位識別方法研究

劉月嬋 ，何元安，商德江, 孫 超

(1.哈爾濱工程大學 水聲技術重點實驗室，哈爾濱 150001；2.船舶系統工程部，北京 100036)

為了突破常規波束形成(conventional beamforming, CBF)在多聲源探測中受陣列物理孔徑的限制所導致的“瑞利限”，眾多高分辨率譜估計法應運而生[1-2]，其中以子空間多重信號分類(MUSIC)算法為代表。該算法通過對基陣接收數據協方差矩陣進行特征分解，從而得到與信號分量相對應的信號子空間與噪聲子空間，利用其正交特性構造出“針狀”空間譜峰，相較于以往的各種方法大大的提高了定位精度及空間分辨率。然而該方法僅能得到聲源方位的估計結果，無法真實反映聲源貢獻的相對大小，不能達到聲源識別的目的。近場聲全息技術[3-6](near field acoustic holography, NAH)作為另一種重要的噪聲源定位識別方法雖然能獲得聲源的真實貢獻量，還可對聲場進行重建進而獲得更豐富的聲場信息，然而該方法在水下輻射噪聲源定位方面的應用較少。究其原因，水下(航行器輻射噪聲通常被視為分布式聲源，且多集中在中低頻段；噪聲源輻射形式未知；噪聲源間相關性未知。這些問題導致以往方法獲得的定位結果缺乏可靠性，并且很難獲得噪聲源輻射強度的精確值。

鑒于特征分解類算法的高分辨率特性，本文提出了一種基于似P范數特征分解的噪聲源定位識別方法。通過特征分解獲得各階特征子空間及子空間響應函數，對各階特征子空間采用預設聲源的方法建立子空間聲源向量重構模型，利用似P范數稀疏性約束[7-9]求解該重構模型，迫使子空間聲源向量能更快，更有效的向真實稀疏源位置收斂，獲得能量更加集中的稀疏解。仿真及水池試驗結果顯示該方法可實現對聲源位置及幅度的準確估計，空間分辨率高。具有很高的工程應用價值。

1 基本理論

(1)

(2)

由于特征子空間與特征矢量具有相同的正交特性，式(2)可寫成：

(3)

式(3)說明陣列對空間聲源的響應可表示成各階特征子空間響應函數的疊加。若Si為第i階特征子空間對應的i階聲源向量，H為聲源向量與接收陣元間的傳遞函數矩陣，那么每階特征子空間均可看成是一個網絡系統，其輸出為子空間響應函數ui：

HSi=ui

(4)

與常規波束形成聲源定位識別方法相同，對于任意階特征子空間，首先將掃描域離散化為N個格點區域，節點處為波束掃描點，各掃描點均存在一預設聲源，該預設聲源可以是點聲源，多極子聲源等，或不同類型聲源的組合形式，文中選定點聲源為預設聲源，其它情況將在后續研究中討論。除目標聲源外，其余為假想聲源，所有預設聲源幅值構成了子空間內N×1維聲源向量Si，認為該子空間響應函數ui是由這些預設源輻射的聲場疊加產生，不同預設聲源類型聯合構成了掃描點與陣元間的M×N維傳遞函數矩陣H。通過求解式(4)可計算出聲源向量Si,即該子空間內波束掃描輸出。考慮式(4)為聲學逆問題，通常采用正則化技術[10-11]獲得問題的近似解：

(5)

其中α稱為正則化參數[12]。其取值與輸入的信噪比有關,通常為矩陣HH-1(或H-1H)最大特征值的0.1%～10%。

由于假想聲源并不真實存在，其幅值在向量Si中表現為零或近似為零，而目標聲源相較與假想源通常具有較大的幅值，這說明Si具有信息局部化的特點，可將其看成是一組稀疏信號。正則化技術雖然解決了逆問題的不適定性，但并不能獲得最優的稀疏解。為了得到能量集中的稀疏解，本文采用基于似P范數的迭代加權算法，將式(4)轉化為如下線性約束優化問題：

(6)

利用Lagrange乘子的方法求解上式,最小化的代價函數為：

(7)

解出關于聲源向量Si和Lagrange乘子λ的形式解如下：

(8)

(9)

綜上得到求解聲源向量Si的迭代公式如下(k≥1)：

(10)

2 仿真分析

為了驗證本文提出的基于似P范數特征值分解的噪聲源定位算法(P-EVD)的有效性,將其聲源定位識別結果與經典的常規波束形成技術(CBF),子空間MUSIC算法及基于幅度補償的MVDR(ac-MVDR)算法[14]進行比較分析。采用由聲壓水聽器構成的7×7均勻網格陣，陣元間距0.3 m，目標聲源位于x-y平面上，掃描域 8 m×8 m，掃描間隔0.25 m，基陣面與掃描面相距z=5 m，設水中聲速為1 500 m/s，背景噪聲為高斯白噪聲，除特別說明，信噪比snr=20 dB，Monte-Carlo實驗數為50，范數p的取值將在后續研究中重點分析，這里不做分析，仿真中取p=1，正則化參數α取最大特征值的1%。

算例1: 兩個非相干單極子聲源A和B分別位于(-1,0,0)和(1,0,0)，發射頻率為1.5 kHz幅值分別為aA=2，aB=1的單頻信號。

圖1 非相干單極子聲源定位

圖1所示是四種算法對非相干單極子聲源的定位結果，由于常規波束形成算法的空間分辨率低，并且旁瓣起伏較大(圖1(a))，從圖中僅能判斷出在(-1,0,0)處存在幅度為1.944的聲源，不能分辨出兩個聲源的空間位置,定位失效；與CBF方法相類似，ac-MVDR算法的空間分辨率雖然有所提高，然而對于文中采用的小孔徑陣列而言，其性能相較于(b)圖的MUSIC算法并無改善；圖1(b)中的MUSIC算法是一種高分辨率空間譜估計算法，其空間分辨率明顯高于CBF，具有尖銳的聚焦峰和低旁瓣級，圖中兩聲源空間位置清晰可見，該方法可實現對非相干聲源的準確定位。但對聲源幅值的估計不盡理想(圖中兩聲源幅值為aA=0.055,aB=0.056)，這是由于MUSIC空間譜估計算法是利用信號子空間與噪聲子空間的正交性估計目標聲源方位，并不能獲得目標聲源的幅值。由上節的闡述可知本文提出的P-EVD算法是以P范數為約束獲得信號能量高度集中的稀疏解，因此在完成對聲源精確定位的基礎上還可實現對聲源幅值的估計，使這一問題得以解決。圖1(d)顯示的是采用P-EVD算法輸出的空間譜，圖中無旁瓣，具有更尖銳的聚焦峰尺度，聲源位置清晰可見，可獲取高精度高分辨率聲源定位結果。圖中顯示兩聲源幅值為aA=1.991，aB=0.991,估計相對誤差為0.67%，估計精度良好，該方法可獲得對聲源方位及幅值滿意的估計結果。

算例2: 考慮兩相干單極子聲源A和B，幅值aA=aB=2，其他聲源參數同算例1。

圖2 相干單極子聲源定位

圖3 分布式聲源定位

從圖2顯示的結果可以看出：① 與算例1的仿真結果相似，CBF的低分辨率使其對兩相干聲源的掃描結果看起來更像是一個聲源，依然無法分辨出聲源的個數及空間位置,定位失效。② 雖然MUSIC算法在非相干聲源定位中展現出高分辨率的優越性，但對于相干聲源，由于無法準確估計聲源個數，因此很難實現對相干聲源的檢測(圖3(b))，且幅度估計結果與實際聲源相差較大，其定位結果同樣失效。③ ac-MVDR具有比CBF和MUSIC算法更低平穩的旁瓣級，然而其空間分辨率不足以區分出兩相距較近的相干聲源。④ P-EVD算法依然具有優異的空間分辨率和精確的幅度估計能力，aA=aB=1.997，其幅度估計相對誤差為0.3%，克服了其他方法不能處理相干聲源的不足。算例1同算例2說明P-EVD方法可實現對聲源位置及幅度的精確估計而無需考慮其相干性。

算例3: 目標聲源為一條有限長線源A和兩個非相干單極子聲源B和C，線源A與x軸夾角為45°，長度1.4 m，其幾何中心位于(-1,1.5,0)，仿真中將其看成是有多個單極子聲源構成的分布式聲源。聲源B與聲源C分別位于 (2,0,0)和(-1,1,0)處，聲源A，B和C均發射頻率為1.5 kHz的單頻信號，幅度為aA=aB=aC=1。聲源布放情況及掃描結果如圖3所示，圖中＊標注處為聲源，□標注處為測量陣列。

算例3仿真結果顯示：① CBF，MUSIC及ac-MVDR算法對分布式聲源的掃描輸出看起來更像是集中聲源，圖中不能準確分辨出聲源數，聲源定位及幅度估計值均存在偏差。MUSIC 與ac-MVDR算法僅僅是消除了CBF中部分旁瓣干擾，依然不能獲得精準的定位結果及幅度估計。② P-EVD的掃描結果與圖中標定的分布式聲源吻合的很好，可同時實現集中式聲源及分布式聲源的定位識別，輸出聲源幅度為aA=0.986,aB=0.979，aC=0.996,估計誤差1.3%。通過該仿真可知，該方法可實現對分布式聲源形態，方向，位置的有效識別及其幅度的準確估計。

算例4：本算例將討論分析頻率及測距陣列孔徑比對本文算法(P-EVD)目標方位及聲源幅度估計的影響，分別給出各自的均方根誤差曲線。頻率變化范圍100～15 000 Hz，測距陣列孔徑比值0.1～30，每個分析頻率及測距孔徑比的Monte-Carlo實驗數為50。

圖4 P-EVD性能分析

圖4給出了P-EVD聲源定位性能受分析頻率及測距孔徑影響的曲線。通過對比分析可以得到如下結論：① P-EVD算法克服了常規波束形成低頻受限的缺點，可實現在低，中，高全頻段內聲源定位，并不受測量距離及陣列孔徑的約束，定位準確有效；② 相對于P-EVD算法對聲源方位估計的寬容性，其聲源幅度估計僅在中高頻段內有效，且測距孔徑比需小于15方可獲得滿意的幅度估計結果。

由于文中采用的是迭代求解計算，且掃描域較大，因此計算速率受到一定程度的制約，表1給出了上述算例中各類算法的計算時間。多次迭代導致P-EVD相對于其他兩種方法計算耗時有所增長，但并不冗長，綜合其在目標方位和幅度估計方面優良性能，該方法依然具有顯著的優勢。

表1 計算時間(單位:秒)

3 水池試驗驗證

為了驗證本文提出算法的可行性，于消聲水池中開展了水下噪聲源定位識別驗證試驗，如圖5所示。由上節對P-EVD算法的性能分析結果，確定試驗中各項參數的選取。試驗中聲源，水聽器陣列及掃描面布放情況如圖6所示。其中掃描面定義為吊放聲源所在的平面，與池底垂直，掃描面大小設定為-4 m≤x≤4 m，-4 m≤y≤4 m，掃描間隔0.25 m。測量面為水聽器陣列所在的平面，與掃描面平行，測量間距z=3 m，考慮到布放條件及水聽器個數的限制，測量系統由10元垂直線陣構成,陣元間距0.2 m，該陣列沿x方向每間隔0.2 m進行一次采樣構成10×10個測點的測量平面。吊放聲源采用球型換能器作為單極子聲源，聲源頻率f=2 kHz，采樣率fs=32.768 kHz，p=1。

圖5 消聲水池及掃描系統

圖6 聲源定位試驗示意圖

圖7 兩非相干單極子聲源試驗結果

工況1：兩非相干單極子聲源A和B入水深度為4 m，x方向間距1 m,其在掃描面上的坐標分布為A：(-0.5,0,0)和B：(0.5,0,0)，聲源B輸出信號幅值是聲源A的2倍，試驗結果以聲源A為基準做歸一化處理，結果如圖7所示。

圖8 兩相干單極子聲源試驗結果

圖7為非相干聲源試驗處理結果。試驗結果顯示P-EVD具有比常規波束形成算法更高的空間分辨率，并且在聲源相距較近時其空間分辨能力要優于子空間MUSIC算法及ac-MVDR算法。圖中所示的聲源幅度估計的相對歸一化結果驗證了本文算法可同時實現目標方位及聲源幅度滿意的識別結果。

工況2:兩相干單極子聲源A和B入水深度為4 m，x方向間距1.6 m,其在掃描面上的坐標分布為A: (-0.8,0,0)和B: (0.8,0,0)，兩聲源輸出信號幅值相同，試驗結果以聲源A為基準做歸一化處理,結果如圖8所示。

圖8中利用似P范數特征值分解(P-EVD)算法對兩相干單極子聲源的試驗結果顯示，與其他方法相比,定位旁瓣得到了明顯的抑制,分辨率高,兩聲源幅值估計結果近似相等. 試驗結果表明P-EVD克服了其他方法不能識別相干聲源的不足，計算過程中可無需考慮聲源間相干性。

4 結 論

提出了一種基于似P范數特征值分解的聲源定位識別方法。該方法通過在掃描點預設不同的聲源類型，認為陣列接受信號由所有預設聲源輻射聲場疊加產生，建立聲源重構模型。利用似P范數獲得最優稀疏解，進而得到高分辨率，背景更為平滑的聲源定位效果。計算機仿真結果及水池試驗結果顯示該方法有效消除了傳統算法分辨率不足，存在柵瓣效應以及相干聲源處理能力不佳的問題，具有良好的背景噪聲抑制能力，可有效實現對單極子聲源及分布式聲源的準確定位及幅度估計而不受聲源間相干性的影響。

[1]Harry L, Van T著.湯俊譯. 最優陣列處理技術[M]. 北京：清華大學出版社，2008.

[2]陳志菲，孫進才，侯宏. 寬帶DOA估計的類MUSIC波束形成算法[J]. 電子學報，2011，39(6)：1257-1260.

CHENG Zhi-fei，SUN Jin-cai，HOU Yong. MUSIC-Like beamforming method for wideband DOA estimation[J]. Acta Electronica Sinica，2011，39(6)：1257-1260.

[3]楊德森，陳歡，時勝國. 基于最大似然估計的輻射噪聲源近場定位方法性能分析[J]. 振動與沖擊，2010，29(12)：67-70.

YANG De-sen，CHEN Huan，SHI Sheng-guo. Orientation measure of radiant noise source at short distance based on maximum likelihood estimation[J]. Journal of Vibration and Shock，2010，29(12)：67-70.

[4]徐亮，畢傳興，陳劍，等. 基于波疊加法的patch近場聲全息及其實驗研究[J]. 物理學報，2007，56(5)：2776-2783.

XU Liang，BI Chuan-xing，CHEN Jian, et al. Algorithm and experimental investigation of patch nearfield acoustic holography based on wave superposition approach[J]. Acta Physica Sinica，2007，56(5)：2776-2783.

[5]張小正，畢傳興，徐亮，等. 基于波疊加法的近場聲全息空間分辨率增強方法[J]. 物理學報，2010，59(8)：5564-5571.

ZHANG Xiao-zheng，BI Chuan-xing，XU Liang, et al. Resolution enhancement of nearfield acoustic holography by the wave superposition approach[J].Acta Physica Sinica, 2010，59(8)：5564-5571.

[6]Christophe L，Manuel M，Alexandre G. Boundary element method for the acoustic characterization of a machine in bounded noisy environment[J]. Journal of the Acoustial Socity of America，2007，121(5)：2750-2757.

[7]Christophe L，Manuel M，Alexandre G. Measurement of confined acoustic sources using near-field acoustic holography[J]. Journal of the Acoustial Socity of America，2009，126(3)：1250-1256.

[8]Figueiredo M A T，Nowak R D. A bound optimization approach to wavelet-based image disconsolation[J]. IEEE Trans Image Processing，2005，2(2)：782-787.

[9]Figueiredo M A T，Bioucas-Dias J M，Nowak R D. Majorization-Minimization algorithms for wavelet-based image restoration[J]. IEEE Trans Image Processing，2007，16(12)：2980-2991.

[10]Chartrand R，Wotao Y. Iteratively reweighted algorithms for compressive sensing[J]. IEEE Trans Acoust，Speech，Signal Processing，2008，34(6)：3869-3872.

[11]Nelson P A，Yooh S H. Estimation of acoustic source strength by inverse methods: Part1 conditioning of the inverse problem[J]. Journal of Sound and Vibration，2000，233(4)：639-664.

[12]Williams E G. Regularization method for near-field acoustic holography[J]. Journal of the Acoustical Society of America. 2001，110(4)：1976-1988.

[13]Candes E J，Wakin M B，Boyd S P. Enhancing sparsity by reweighted l1 minimization[J]. Journal of Fourier Analysis and Application. 2008，14(5-6)：877-905.

[14]陳歡，何良，楊德森，等.基于幅度補償的MVDR水下噪聲源近場定位識別方法研究[J]. 振動與沖擊，2012，31(2)：51-54.

CHEN Huan，HE Liang，YANG De-sen，et al. Underwater noise sources identification in near-field location based on MVDR method with amplitude compensation [J]. Journal of Vibration and Shock，2012，31(2)：51-54.