基于不動點迭代的增廣拉格朗日聲源識別算法

賀巖松， 賈晨陽， 黃琳森， 昝 鳴， 徐中明

(1.重慶大學(xué) 汽車工程學(xué)院，重慶 400030； 2. 重慶大學(xué) 機(jī)械傳動國家重點實驗室，重慶 400030)

近場聲全息是聲源識別和聲學(xué)成像的重要方法之一，其發(fā)展和應(yīng)用十分廣泛[1-2]。該方法可以在近場記錄得到高空間頻率的倏逝波，從而突破瑞利判據(jù)的限制，提高聲場重建性能。

早在20世紀(jì)80年代，Williams等[3-4]就首次提出基于空間傅里葉變換的近場聲全息理論，但該方法對聲源形狀和全息面有諸多限制，且進(jìn)行傅里葉變換時會產(chǎn)生卷繞誤差。為突破該方法的限制，先前學(xué)者提出了許多新的方法，如邊界元法[5-6]、等效源法[7-9]、統(tǒng)計最優(yōu)[10-11]等。其中等效源法可以應(yīng)用于任意形狀聲源，且擁有魯棒性好、效率高、易于實現(xiàn)等諸多優(yōu)點。等效源法主要思想是將實際聲源等效為許多虛擬點源的疊加，通過建立聲源面到全息面的傳遞關(guān)系，求解得到源強矢量，從而重建聲場。在應(yīng)用等效源法求解逆問題時，會產(chǎn)生不適定性問題，需要對重建過程進(jìn)行正則化處理。傳統(tǒng)方法主要應(yīng)用L2范數(shù)正則化方法，包括Tikhonov正則化[12-13]和截斷奇異值分解[14-15]等，該方法主要針對低頻聲源，對中高頻聲源聲場重建效果較差。

近年來，基于L1范數(shù)最小化問題的稀疏正則化方法被逐漸引用到聲源識別中。Chardon等[16]首次將壓縮感知理論引入近場聲全息，并構(gòu)建近場聲全息的稀疏框架。Gade等[17-18]提出寬帶聲全息(wideband acoustic holography, WBH)，在中高頻擁有很好的重建性能，在一定程度上彌補了傳統(tǒng)方法的短板，但在低頻時對相干聲源的識別效果不佳。Fernandez-Grande等[19]在稀疏等效源法的框架下應(yīng)用CVX工具箱進(jìn)行求解，在較寬的頻率范圍內(nèi)擁有較好的重建精度，但該方法比較耗時，效率不高。張晉源等[20]在傳統(tǒng)壓縮波束形成的基礎(chǔ)上，應(yīng)用迭代重加權(quán)L1范數(shù)進(jìn)行改進(jìn)，從而有效地提高了聲源的空間分辨率和旁瓣衰減能力。Hald[21]闡述并對比五種不同的迭代稀疏等效源法，包括寬帶聲全息和快速迭代收縮閾值算法(fast iterative shrinking threshold algorithm, FISTA)，并在仿真和試驗的基礎(chǔ)上，討論了不同方法的優(yōu)化選擇。樊小鵬等[22]提出分裂增廣拉格朗日反卷積聲源識別算法，利用交替迭代乘子法求解反卷積問題，顯著提高了聲源識別精度和計算效率。增廣拉格朗日方法(augmented Lagrangian method, ALM)是解決L1范數(shù)正則化問題的重要方法，該方法計算效率高，且對懲罰參數(shù)有很好的魯棒性，在壓縮感知領(lǐng)域有諸多應(yīng)用[23-25]。

本文致力于在等效源法的稀疏框架之下，為提高聲源識別的頻率范圍和分辨率，提出一種基于增廣拉格朗日方法的聲場重建算法；為證明該算法在全頻帶下的重建性能，仿真時，將該算法與Tikhonov正則化方法，寬帶聲全息和快速迭代收縮閾值三種算法進(jìn)行對比，同時檢驗其對不同信噪比(signal-to-noise ratio, SNR)和全息距離的適應(yīng)性。最后，通過試驗驗證該算法的性能。

1 基本理論

1.1 等效源法

等效源法的基本原理如圖1所示，其將振動體(聲源)產(chǎn)生的聲場由一系列等效點源產(chǎn)生的聲場疊加代替，而這些等效源的強度可以利用全息面的聲壓反求得到，從而實現(xiàn)聲場重建和聲源的準(zhǔn)確定位。

如圖1所示，假設(shè)全息面有M個測量點，等效源面分布有N個虛擬等效點源，則第m個測量點的聲壓可表示為

(1)

式中：qn為第n個等效源的源強;g(rm,rn)為第n個等效源和第m個測量點之間的聲壓傳遞函數(shù)，表達(dá)式為

(2)

圖1 等效源法示意圖Fig.1 Schematic diagram of equivalent source method

p=Aq

(3)

式中，q的表達(dá)式為q=iρckQ，Q為體積速度向量。通過對式(3)逆向求解可得到等效源強度，并可以重建得到重建面聲壓，即

pR=Gq

(4)

式中：pR為重建面聲壓;G為等效源面到重建面的傳遞矩陣。

考慮到測量成本等因素，實際應(yīng)用中全息面?zhèn)髀暺鞯臄?shù)量遠(yuǎn)遠(yuǎn)小于等效源的數(shù)量，即M?N，因此，在進(jìn)行欠定方程的逆向求解時會產(chǎn)生不適定問題。為解決此問題，Tikhonov正則化方法是最為常用的方法之一，其引入L2范數(shù)對源強向量進(jìn)行約束，可表示為

(5)

式中，λ為正則化參數(shù)，對其合理選擇可以很好地平衡量化誤差和源強向量。

當(dāng)?shù)刃г吹姆植际窍∈杌蚪葡∈钑r，對源強向量進(jìn)行稀疏約束能夠更加精確地得出聲源信息。理想情況下，可以應(yīng)用L0范數(shù)約束求得等效源強的稀疏解，此時式(5)轉(zhuǎn)化為

(6)

(7)

1.2 ALM-FP算法

在壓縮感知理論框架下，L1范數(shù)正則化方法成為稀疏重建問題的重要方法。本文應(yīng)用如下L1范數(shù)正則化模型求解稀疏重建逆問題，在選擇合理的正則化參數(shù)的情況下，該問題與式(7)等效

(8)

增廣拉格朗日方法對懲罰參數(shù)有很好的魯棒性，且易于求解，同時應(yīng)用不動點(fixed point, FP)迭代方法求解其子問題，使該算法擁有很好的全局收斂性，因此本文應(yīng)用該方法解決上述最小化問題。

增廣拉格朗日方法求解過程為

(9)

(10)

對于式(10)的拉格朗日最小化問題，我們應(yīng)用不動點迭代進(jìn)行求解，令

(11)

則式(10)轉(zhuǎn)化為

(12)

根據(jù)文獻(xiàn)[26]，該問題求得的結(jié)果為

式中:τ為閾值因子，其取值范圍為0<τ≤2/λmax，λmax為矩陣ATA的最大特征值，本文令τ=1.8/λmax；?g(q)的表達(dá)式為

(14)

以上即為ALM-FP算法實現(xiàn)的整個過程，其具體的算法步驟以偽代碼形式總結(jié)如下：

計算全息面聲壓p以及傳遞矩陣A和G

初始化，令λ0=0，q0=0，k=0，閾值因子τ=1.8/λmax。

通過式(14)求得?g(qk)；

通過式(13)計算qk+1；

通過式(9)更新拉格朗日算子λ；

更新迭代次數(shù)k=k+1；

把求解得到的q當(dāng)作下一步迭代的初始值

End

輸出迭代完成得到的最優(yōu)q值；

重建得到重建面聲壓;pR=Gq。

2 仿真分析

為驗證ALM-FP算法的聲場重建性能，對該算法進(jìn)行相應(yīng)的仿真分析。仿真時聲源、重建面和全息面按照圖2(a)布置。定義聲速為340 m/s，空氣密度為1.29 kg/m3，目標(biāo)聲源采用脈動球源，其半徑設(shè)為0.01 m，振速為2.5×10-2m/s；設(shè)定等效源面和重建面到聲源面距離分別為0.001 m和0.02 m，全息距離為0.3 m，等效源面分布有41×41個等效源，網(wǎng)格間距為0.02 m，重建面網(wǎng)格尺寸為0.8 m×0.8 m，網(wǎng)格間距為0.02 m。為保證仿真和試驗的一致性，仿真采用36通道圓形隨機(jī)陣列，陣列點的坐標(biāo)如圖2(b)所示。所有的仿真測試均在Windows10操作系統(tǒng)的電腦(AMD A10-5800K 處理器，3.80 GHz，8.00 GB RAM)上的MATLAB R2016a(64位)中完成。

圖2 仿真布置圖和傳聲器陣列Fig.2 Layout diagram of simulation and microphone array

本次仿真分析以單聲源和相干雙聲源兩種狀況進(jìn)行分析，以Tikhonov正則化、WBH、FISTA作為新算法的對比研究對象。單聲源時，聲源位于坐標(biāo)(0,0,0)m處，即為聲源表面的正中心，雙聲源時，聲源的位置分別為(0.2, 0, 0)m和(-0.2, 0, 0)m。另外，為了定量地分析聲場重建精度，定義幅值重建誤差E1和平均相位誤差E2如下

(15)

(16)

式中：pcal為計算得到的重建面聲壓；pt為重建面理論聲壓;L為重建面網(wǎng)格點數(shù)；pcal,i和pt,i分別為重建面第i個網(wǎng)格點的計算聲壓和理論理論聲壓；φ(·)為取相位。幅值重建誤差E1和平均相位誤差E2越小，表明重建精度越高。

2.1 性能表現(xiàn)

為檢驗ALM-FP算法的性能，在單聲源和相干聲源情況下對該算法進(jìn)行仿真模擬，并與Tikhonov正則化、WBH、FISTA三種方法進(jìn)行對比。仿真時加入信噪比為30 dB的高斯白噪聲，動態(tài)顯示范圍為15 dB。

表1和表2分別展示了四種方法在500 Hz，1 500 Hz以及3 000 Hz三個頻率下二維重建面聲壓云圖，其中表1為單聲源，表2為相干聲源。從表中可以看出，在兩種識別過程中， Tikhonov正則化方法在三個頻率下的主瓣寬度都比較大，在低頻時尤為明顯，真實聲源無法被準(zhǔn)確識別和定位，這主要是由于測點數(shù)目較少導(dǎo)致；WBH方法對單聲源識別表現(xiàn)良好，在低頻時識別相干聲源失效，無法識別出兩個聲源位置。FISTA和ALM-FP兩種方法都能得到很好的結(jié)果，但相對而言，ALM-FP方法在識別過程中擁有更好的分辨率，而且其最大聲壓值與理論最大聲壓最為接近。

表1 單聲源重建面聲壓云圖(仿真)Tab.1 Sound pressure maps of reconstruction surface for single sound source(simulation)

表2 相干聲源重建面聲壓云圖(仿真)Tab.2 Sound pressure maps of reconstruction surface for coherent sound source(simulation)

三個頻率下四種方法的重建面中間行聲壓對比圖，如圖3所示。圖3(a)～圖3(c)為單聲源，圖3(d)～圖3(f)為相干聲源。其結(jié)果基本印證了上述結(jié)論，三個頻率下，Tikhonov正則化方法得到的聲壓值都與理論值相差較大；單聲源時，ALM-FP方法與理論值契合度最高，效果最好；相干聲源狀況下， WBH在低頻時基本失效，F(xiàn)ISTA和ALM-FP表現(xiàn)良好，同樣是ALM-FP方法與理論值更為接近。

圖3 重建面中間行聲壓Fig.3 Sound pressure in the middle row of the reconstruction surface

表3為相干聲源情況下四種算法在1 500 Hz時的計算時間，該結(jié)果為十次計算下的平均值，其中FISTA和ALM-FP的迭代次數(shù)均為1 000次。從結(jié)果可以看出，Tikhonov和WBH算法擁有較好的計算效率，其中WBH計算時間在0.1 s以下，計算效率最高，但兩種算法在精度上都有一定的缺陷。FISTA和ALM-FP算法的計算時間分別為21.077 s和10.359 s，相對其他兩種算法雖然效率不高，但仍在可接受的范圍內(nèi)，而且，相對FISTA算法，ALM-FP在精度和計算效率上都具有一定的優(yōu)勢。

表3 四種方法計算時間(相干聲源1 500 Hz)Tab.3 Computational efficiency comparison of the four methods (coherent sound source at 1 500 Hz)

2.2 誤差分析

通過定量地分析算法的幅值誤差和相位誤差，能夠更好的看出算法在全頻段下的聲場重建性能，改變仿真時的信噪比和全息距離，能夠直觀的看出各算法對信噪比以及測量距離的適應(yīng)性。

圖4分別展示了單聲源和相干聲源兩種情況下，從200～4 000 Hz四種算法的重建面幅值誤差曲線。每條曲線每隔50 Hz取一個頻率點，單個頻率的誤差計算10次以取得平均值。單聲源時，Tikhonov正則化算法在低頻時的幅值重建精度要優(yōu)于中高頻，但整體上重建誤差在35%以上，相對其他算法表現(xiàn)較差。WBH、FISTA和ALM-FP三種算法均有比較好的幅值重建精度，重建誤差均在20%以下，且中高頻的重建效果要好于低頻。ALM-FP算法整體上的重建誤差在10%以下，相對其他算法在低頻的重建誤差更小，對低頻單聲源聲場重建更有優(yōu)勢。雙聲源時，Tikhonov正則化算法與單聲源時表現(xiàn)類似。WBH算法的曲線上存在一個明顯的轉(zhuǎn)換頻率，在轉(zhuǎn)換頻率以下的低頻時出現(xiàn)較大的重建誤差，大于轉(zhuǎn)換頻率時，重建誤差較小。FISTA和ALM-FP依然擁有較高的幅值重建精度，相對WBH和Tikhonov正則化優(yōu)勢明顯，但整體上而言，特別是在低頻時，ALM-FP算法擁有更小的重建誤差，在進(jìn)行聲場重建時效果更好。

圖4 聲壓幅值重建誤差曲線Fig.4 Curves of source pressure amplitude reconstruction error

為了進(jìn)一步分析算法的重建性能，以相干聲源為例，圖5給出了四種算法在1 500 Hz時的重建面網(wǎng)格點的相位以及四種算法在200～4 000 Hz頻率段的平均相位誤差，由于重建面網(wǎng)格點較多，隨機(jī)選取連續(xù)的100個點為代表以體現(xiàn)重建聲壓的相位情況?？梢园l(fā)現(xiàn)，Tikhonov正則化算法重建得到的聲壓相位與理論相位有較大的偏移，相位的重建效果較差，其他三種方法得到的結(jié)果與理論結(jié)果有相同的走勢，其中WBH有些許的偏移，F(xiàn)ISTA與ALM-FP能夠與理論曲線很好的契合，說明其能夠得到很好的相位結(jié)果。通過平均相位誤差曲線可以看出，與相干聲源幅值誤差曲線有相似之處，Tikhonov正則化算法會產(chǎn)生較大的相位誤差，特別是中高頻段。WBH低頻相位誤差較大，高于轉(zhuǎn)換頻率時效果較好。FISTA與ALM-FP在整個頻段都有較小的相位重建誤差，F(xiàn)ISTA算法誤差略小一點，但整體而言，兩種算法均表現(xiàn)出很好的相位重建性能。

圖5 相位分布及相位誤差Fig.5 Phase distribution and phase error

通過算法在單雙聲源的重建面幅值誤差曲線可以發(fā)現(xiàn)，不同頻率時，當(dāng)單聲源重建較好時，雙聲源的重建效果并不好，在相干雙聲源情況下，對算法的要求更苛刻，因而該情況下，也更能反映算法對信噪比和全息距離的適應(yīng)性。四種方法在不同信噪比的情況下的重建誤差云圖，如圖6所示，信噪比范圍為10～60 dB，全息距離為0.3 m。四種方法在不同全息距離下的重建誤差云圖，如圖7所示，全息距離范圍為0.05～0.50 m，信噪比為30 dB。兩圖均在相干聲源情況下計算所得，圖中顏色條的范圍為0～100 %。

圖6 不同信噪比下各方法的重建誤差云圖Fig.6 Reconstruction error map of methods with different SNR

圖7 不同全息距離下各方法的重建誤差云圖Fig.7 Reconstruction error map of methods with different hologram distance

從圖6可以看出，Tikhonov正則化方法在低頻重建效果優(yōu)于中高頻，且中高頻時基本失效，低頻時，信噪比較高時的重建誤差較小。WBH方法在低頻失效，在中高頻對信噪比的適應(yīng)性較好，在不同信噪比下都有比較好的重建效果。FISTA與ALM-FP兩者有相似的趨勢，但ALM的整體效果更好，其重建誤差隨頻率升高和信噪比增大而減小，只在頻率小于800 Hz、信噪比小于15 dB時，重建誤差較大。

從圖7可以看出，全息距離的減小有利于聲場重建的效果。圖7(a)為Tikhonov正則化方法，通過減小全息距離，其重建誤差可以減小為原來的一半。另外減小全息距離可以拓寬WBH方法的頻率范圍，在全息距離較小時，WBH在低頻時也可以有較好的重建精度。ALM-FP方法對全息距離有較好的適應(yīng)性，在不同全息距離下，都有比較好的重建精度，只在全息距離接近0.5 m的低頻時，有較大的重建誤差。FISTA方法整體上在不同全息距離下有比較好的重建精度，但頻率適用性上要劣于ALM-FP方法。

3 試驗驗證

驗證算法性能的試驗布局圖如圖8所示，將單頻信號激勵的音響作為聲源，采用B & K公司、直徑0.65 m的36通道扇形輪陣列測量聲壓信號，陣列集成4958型傳聲器，各傳聲器測量的聲壓信號經(jīng)PULSE 3660D型數(shù)據(jù)采集系統(tǒng)同時采集并傳輸?shù)絇ULSE LABSHOP中進(jìn)行數(shù)據(jù)分析，信號采樣頻率為32 768 Hz，將得到的數(shù)據(jù)在MATLAB軟件中進(jìn)行計算，得到重建面聲壓云圖。圖8(a)為單聲源布局，聲源坐標(biāo)為(0, 0, 0)m，圖8(b)為相干聲源布局，聲源坐標(biāo)為(0.2, 0, 0)m、(-0.2, 0, 0)m，聲源與陣列之間的距離均為0.2 m，其他參數(shù)設(shè)置與仿真時相同。

圖8 試驗布局圖Fig.8 Layout of experiment

與仿真相對應(yīng)，表4和表5分別展示了試驗時單聲源和相干雙聲源在500 Hz，1 500 Hz和3 000 Hz三個頻率下的重建面聲壓云圖，其結(jié)果能很好地驗證仿真的結(jié)論。單聲源時，但是Tikhonov正則化方法雖然能識別到噪聲源的位置，但是出現(xiàn)較多的鬼影，而且主瓣較大，分辨率很低。WBH、FISTA和ALM-FP三種方法能比較好識別噪聲源的位置，且隨著頻率的升高，聲源主瓣變小，分辨率變好，但整體上特別是在低頻條件下，ALM-FP較其他兩種方法有更好的分辨率，與仿真的結(jié)果相對應(yīng)。相干雙聲源情況下，Tikhonov正則化方法出現(xiàn)較多的旁瓣和鬼影，基本無法判別噪聲源的位置，聲源識別效果很差。WBH方法在低頻時失效，將兩個聲源識別成為中間的一個聲源，在中高頻能很好地識別噪聲源的位置。FISTA和ALM-FP兩種方法能很好地識別噪聲源的位置，不過會出現(xiàn)一些微小的旁瓣，原因可能是由于試驗是在普通室內(nèi)進(jìn)行的，壁面的聲波反射無法避免，從云圖的結(jié)果來看，ALM-FP方法在三個頻率下的分辨率更高，而且能夠更好地消減旁瓣，效果更好。

表4 單聲源重建面聲壓云圖(試驗)Tab.4 Sound pressure maps of reconstruction surface for single sound source(experiment)

表5 相干聲源重建面聲壓云圖(試驗)Tab.5 Sound pressure maps of reconstruction surface for coherent sound source(experiment)

4 結(jié) 論

基于等效源法和壓縮感知理論，本文提出一種基于增廣拉格朗日方法的聲源識別算法ALM-FP，并給出了該方法的理論推導(dǎo)，模擬仿真了單聲源和相干雙聲源情況下的聲源識別成像效果和隨頻率變化的誤差曲線，給出了不同信噪比和全息距離下的重建誤差云圖，最后通過試驗驗證了仿真結(jié)果。

通過仿真和試驗結(jié)果表明，ALM-FP算法能夠在較寬的頻率范圍內(nèi)進(jìn)行準(zhǔn)確的聲源識別和聲場重建，對比其他三種方法，該方法擁有更好的重建精度，聲源識別分辨率更高，特別是在低頻區(qū)域，其優(yōu)勢更加明顯。另外該方法對信噪比和全息距離擁有很好的適應(yīng)性，但是在測量環(huán)境特別惡劣時，比如低頻時信噪比很小或全息距離很大，會嚴(yán)重影響該方法的重建效果，產(chǎn)生較大的重建誤差。