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

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

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

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

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

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

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

1 基本理論

1.1 等效源法

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

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

(1)

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

(2)

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

p=Aq

(3)

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

pR=Gq

(4)

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

考慮到測量成本等因素，實際應用中全息面傳聲器的數量遠遠小于等效源的數量，即M?N，因此，在進行欠定方程的逆向求解時會產生不適定問題。為解決此問題，Tikhonov正則化方法是最為常用的方法之一，其引入L2范數對源強向量進行約束，可表示為

(5)

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

當等效源的分布是稀疏或近似稀疏時，對源強向量進行稀疏約束能夠更加精確地得出聲源信息。理想情況下，可以應用L0范數約束求得等效源強的稀疏解，此時式(5)轉化為

(6)

(7)

1.2 ALM-FP算法

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

(8)

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

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

(9)

(10)

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

(11)

則式(10)轉化為

(12)

根據文獻[26]，該問題求得的結果為

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

(14)

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

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

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

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

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

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

更新迭代次數k=k+1；

把求解得到的q當作下一步迭代的初始值

End

輸出迭代完成得到的最優q值；

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

2 仿真分析

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

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

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

(15)

(16)

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

2.1 性能表現

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

表1和表2分別展示了四種方法在500 Hz，1 500 Hz以及3 000 Hz三個頻率下二維重建面聲壓云圖，其中表1為單聲源，表2為相干聲源。從表中可以看出，在兩種識別過程中， Tikhonov正則化方法在三個頻率下的主瓣寬度都比較大，在低頻時尤為明顯，真實聲源無法被準確識別和定位，這主要是由于測點數目較少導致；WBH方法對單聲源識別表現良好，在低頻時識別相干聲源失效，無法識別出兩個聲源位置。FISTA和ALM-FP兩種方法都能得到很好的結果，但相對而言，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)為相干聲源。其結果基本印證了上述結論，三個頻率下，Tikhonov正則化方法得到的聲壓值都與理論值相差較大；單聲源時，ALM-FP方法與理論值契合度最高，效果最好；相干聲源狀況下， WBH在低頻時基本失效，FISTA和ALM-FP表現良好，同樣是ALM-FP方法與理論值更為接近。

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

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

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

2.2 誤差分析

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

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

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

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

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

通過算法在單雙聲源的重建面幅值誤差曲線可以發現，不同頻率時，當單聲源重建較好時，雙聲源的重建效果并不好，在相干雙聲源情況下，對算法的要求更苛刻，因而該情況下，也更能反映算法對信噪比和全息距離的適應性。四種方法在不同信噪比的情況下的重建誤差云圖，如圖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正則化方法在低頻重建效果優于中高頻，且中高頻時基本失效，低頻時，信噪比較高時的重建誤差較小。WBH方法在低頻失效，在中高頻對信噪比的適應性較好，在不同信噪比下都有比較好的重建效果。FISTA與ALM-FP兩者有相似的趨勢，但ALM的整體效果更好，其重建誤差隨頻率升高和信噪比增大而減小，只在頻率小于800 Hz、信噪比小于15 dB時，重建誤差較大。

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

3 試驗驗證

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

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

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

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

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