基于L-曲線參數優化的均勻聲場重建算法

張 旭 趙洪亮 程希光 裴景春

(1 山東科技大學電子信息工程學院 青島 266590)

(2 山東科技大學電氣與自動化工程學院 青島 266590)

0 引言

揚聲器陣列的均勻聲場重建是指控制揚聲器陣列輸入的幅值和相位，使得控制區域的重建聲場均勻化，同時增大明區能量占比。為達到目的，已提出多種控制方法，大體可分為兩類：一類將建立能量控制作為目標函數，如能量對比度控制(Acoustic contrast control,ACC)法[1?2]、聲能量差最大化控制(Acoustic energy difference maximization,AEDM)法[3]等；另一類將聲場重建誤差作為目標函數，如基于最小二乘準則的多點匹配法(Pressure matching with the least-squares criterion,PM-LS)[4?5]等。

在聲場重建中，往往出現不適定問題[6]。其中最常用的是使用Tikhonov 正則化來解決不適定的聲場問題[7]，而正則化程度取決于正則化參數。對于正則化參數的選取，文獻[8–9]中提出了使用平滑矩陣的廣義Tikhonov 正則化。在文獻[8]中，導數平滑范數應用于自適應波場合成，以減少揚聲器激勵信號的相位變化。Betlehem 等[2]設計了一種方法，計算揚聲器權重的平方和來限定Tikhonov參數上限，然后采用牛頓法確定正則化參數。另一個方法就是截斷奇異值分解(Truncate singular value decomposition,TSVD)，對傳遞矩陣進行奇異值分解，并且只有較重要的特征值才會被用于聲場重建。該方法減少并消除了較小奇異值對解的影響，但是對于正則化參數的選取很難確定一個合理的標準[10]。同時揚聲器陣列功率受到物理系統的制約，需要考慮揚聲器驅動信號功率與重現誤差之間的平衡問題。因此怎樣在考慮揚聲器功率限制的同時合理地選擇正則化參數是亟需解決的問題。

在本文中，聲場重建模型以最小化聲場重建誤差為目標函數，并約束揚聲器功率上限。針對正則化參數選擇問題，將L-曲線法[11]引入均勻聲場重建，該方法以重建誤差作為橫軸，揚聲器功率作為縱軸得到擬合曲線，然后選取該曲線上曲率最大的點所對應的參數值作為Tikhonov 正則化參數的選值。作為對比，分別用傳統最小二乘法、廣義交叉驗證(Generalized cross-validation,GCV)法以及L-曲線法進行性能仿真，并在限定功率下對明區重建聲場進行測試。仿真及實驗證明，L-曲線法平衡了揚聲器功率與重建誤差，且陣列能量效率較高，具有良好的聲場重建效果。

1 基于L-曲線法的均勻聲場重建算法

1.1 均勻聲場重建模型

將控制區域劃分為明區以及暗區。明區為目標聲場，為達到良好的均勻聲場，該期望聲壓幅值設為1。基于最大化控制目標區域聲場能量的目的，設置暗區并將暗區期望聲壓設為0。聲場重建區域明暗區劃分如圖1所示。

圖1 聲場重建區域明暗區示意圖Fig.1 Schematic diagram of the light and dark area of the sound field reconstruction

在控制區域0～B3內，明區范圍B1～B2設置M個控制點，暗區0～B1以及B2～B3設置N個控制點。明暗區通過控制點到陣元的距離(m=1,2,··· ,Mb)和(m=1,2,··· ,Nd)得到相應的陣列響應和而對于L個揚聲器陣元權重gL，兩者關系表示如下：

式(1)中，Hb為明區傳遞矩陣，Hd為暗區傳遞矩陣。傳遞矩陣可定義為

整個控制區的傳遞矩陣Hc=[Hb,Hd]T，控制區期望聲壓Pc=[Pb,Pd]T。

本文優化模型是在考慮最小化控制區域重建聲場的誤差函數的基礎上，約束揚聲器權重的平方和上限，其優化模型如下：

其中，∥·∥為2-范數，g為揚聲器權向量，K1為其揚聲器功率限制參數。

1.2 L-曲線選取正則化參數

為獲得揚聲器陣列權重，需對優化模型進行求解。式(3)具有二次目標和單個二次約束，因此該優化模型為凸優化問題。式(3)對應的廣義拉格朗日罰函數可以寫為

式(4)中，λ為揚聲器陣列約束權重因子，用于平衡重建聲場誤差與揚聲器權重g，λ0。式(4)可通過求解，得到揚聲器權向量如下：

式(5)中，H 為共軛轉置，I為單位矩陣。矩陣Hc不是滿秩矩陣，因此進行奇異值分解(Singular value decomposition,SVD)，如下：

其中，Σ ∈Rn×n為對角矩陣，Σ=diag(s1···sn),s1s2···sn0。U ∈Rm×n和V ∈Rn×n分別為m×n和n×n階酉矩陣：

因此，式(5)可以改寫為

以μ=lg∥Hcg(λ)?Pc∥2為橫軸、v=lg∥g(λ)∥2為縱軸并以λ為參量得到擬合曲線，該曲線上曲率最大的點所對應的參數，就是L-曲線法選擇的正則化參數λL，該曲率最大的位置即為L-Corner，其表達式為

其中，μ′、v′和μ′′、v′′分別是μ、v對λ的一階、二階導數。將求得的參數值λL帶入式(5)得到揚聲器陣列權向量gL。

1.3 聲場均勻控制的性能指標

為衡量揚聲器陣列在控制區域的均勻聲場重建效果，本文從以下三個參數對重建聲場性能進行評價，分別為明區重建均方誤差、陣列能量效率以及聲場均勻度。重建性能指標如下。

1.3.1 明區重建均方誤差

定義為明區重建聲場與期望聲場的聲壓幅值之差的均方和與期望聲場的聲壓幅值均方和之比，如式(13)所示：

式(13)中，pb、Pb分別為明區重建聲場聲壓與明區期望聲壓；Hb為明區傳遞函數矩陣。該指標用于衡量聲場重建效果的準確度。

1.3.2 聲能量陣列效率

除了明區重建誤差外，仿真中還需要比較各種情況下揚聲器陣列輻射到目標區域的聲能量大小，本文使用揚聲器陣列在明區的聲能量效率進行評價，明區控制點上平均聲能量可以表示為

式(14)中，Q為目標區域控制點的數量，ptm為明區中第m個控制點的陣列響應，Hb為對應的揚聲器陣列單元到明區內控制點的傳遞函數矩陣。陣列能量效率的定義為明區M個控制點的平均聲能量與揚聲器陣列輸出總能量的平均聲能量的比值，其表達式如下：

該方程代表了揚聲器陣列輻射到明區的聲能量效率，即對總能量的利用效率。方程中對控制區域控制點的能量取均值是為了避免控制點數量的變化對該參數的影響。

1.3.3 明區聲壓級均勻長度和聲壓級波動標準差

對聲場均勻度進行評價時，主要采用明區的均勻長度和聲壓級波動的標準差兩個參數。明區的均勻長度定義為明區聲壓幅值與均值相差少于變動范圍3 dB 內的距離長度，該指標可以反映出明區中均勻聲場所占的比例。而明區聲壓級波動的標準差則可以反映出重建區域整體的聲場波動情況。明區的均勻長度rT和聲壓級波動的標準差Sd分別表示為

其中，表示為明區聲壓向量的均值；pb(rT)為符合條件的明區聲壓向量；M為明區的控制點數量。

2 數值仿真

2.1 L-曲線法的參數仿真

如圖1 所示，取揚聲器陣元個數L=16，激勵頻率f=1000 Hz，陣元間距d=0.1 m，陣列高度H1=7 m。在控制區0～60 m 內，設置10～50 m為明區范圍，即保持聲場均勻的區域，其余為暗區范圍。控制區期望聲場Pc=[Pb,Pd]H。Pb、Pd分別為明區期望聲壓和暗區期望聲壓。基于暗區聲壓數值盡量小的考慮，因此設置Pd=[0,0,··· ,0N]H。明區期望聲壓表示如下：

其中，r=[r1,r2,...,rM]為陣元中心到明區控制點的距離。對重建模型根據L-曲線法求解λ擬合圖，如圖2所示。

由圖2 可以看出f=1000 Hz 時，擬合曲線有一個非常明顯的拐角(L-Corner)，該拐角位于L-曲線的垂直部分與水平部分相交的位置。水平分布所對應的正則化解主要由明區重建誤差主導，而垂直部分由揚聲器權重主導。因此，L-曲線的這個拐角為揚聲器權重與明區重建誤差都較小的一個平衡點。因而在拐角處得到參數值，即λL=0.37874。相應地，在陣列頻率為0～4000 Hz 時，λL對應的取值如圖3所示。可以看出，λL在低頻段取值在0.4左右，在頻率為3500 Hz 以上的高頻段最優取值為1左右。

圖2 L-曲率擬合圖Fig.2 L-curvature fit map

圖3 λL 取值Fig.3 λLValue

2.2 重建性能指標仿真

激勵頻率范圍設置為0～2000 Hz，為有效重現該范圍內的聲場，選取陣元間距為d=0.1 m。根據定義的性能指標，在0～2000 Hz 頻帶上，分別對傳統最小二乘法、GCV 法以及L-曲線法進行仿真比較。

如圖4(a)、圖4(b)所示，可以明顯看出在低頻段約0～600 Hz，未正則化的最小二乘法聲場重建均勻度存在明顯的波動，而高頻段聲場均勻度隨著頻率的變化較為連續且對應的曲線較為平滑。陣列能量效率在0～600 Hz 頻段內較低，波動范圍為?160～?140 dB，而在600～1400 Hz 間，陣列效率頻率上升而提升至?20 dB。基于以上分析可以看出，傳統的最小二乘法存在低頻段能量過低以及均勻度波動問題。

圖4 重建性能指標仿真對比結果Fig.4 Reconstruction performance index simulation comparison result

為了改善該情況，需要對其進行正則化處理。對于正則化因子的選擇，仿真選取了GCV 法與L-曲線法進行對比。GCV 法選擇最佳正則化參數就是選取使GCV 函數Gλ最小時對應的λ值。從圖4中可以看出，基于GCV 法選擇的正則化參數平滑了低頻段的波動，均勻度穩定且隨著頻率增加而上升，重建誤差在20%以下。然而圖4(d)顯示GCV法由于追求重建誤差，陣列能量效率沒有顯著提升。

GCV 法與L-曲線法皆可以對傳統的最小二乘法低頻段波動進行優化，相比于GCV 法，L-曲線法的陣列能量效率整體提升了130 dB,極大改善了低頻段的陣列能量效率過低的問題。均勻長度方面如圖4(a)、圖4(b)所示，L-曲線法在0～1000 Hz 頻段與1400～2000 Hz 頻段分別保持一致，并且整個頻段聲壓級標準差變化平緩，波動在0.2～0.25 之間。基于L-曲線法的明區重建誤差在25%，整體均勻只在頻率300 Hz處有一個波動。圖4(c)結果表示基于L-曲線法的重建誤差低于GCV 法的重建誤差5%，該原因是L-曲線法約束了揚聲器功率，犧牲了一定程度的重建誤差。GCV 法為了到達良好的重建誤差，增加了揚聲器陣列輸出，因此GCV 法陣列能量效率較低。在實際系統實現中，揚聲器功率是極為重要的指標，因此L-曲線法平衡了揚聲器功率與重建誤差，更為符合實際要求。

3 實驗比較

本節將揚聲器陣列未控制的聲場效果與GCV法及L-曲線法的重建聲場進行比較。本實驗將明區聲場縮小為1～5 m，陣元中心距離地面為2 m，揚聲器陣列數量為16個，陣元間距為0.6 m，測量間距為0.2 m，如圖5所示。期望聲壓設為85 dB，并將三種方法的揚聲器功率均設置為以L-曲線功率為標準進行計算。實驗以及仿真頻率均取1000 Hz，在此頻率下，L-曲線法取正則化參數值為λL=0.4480，GCV法取參數值為λG=2.66×10?5。

圖5 實驗裝置圖Fig.5 Experimental setup diagram

首先對L-曲線法及GCV 法進行均勻聲場重建仿真。如圖6所示，該聲場左側聲壓最大值處為揚聲器陣列，陣列根據求得的揚聲器權重偏轉一定角度，在明區范圍內進行均勻聲場重建。由圖6 的聲壓級分布圖可以看出，L-曲線法在明區重建區域誤差以及均勻度略差于GCV 法。然而，GCV 法的陣列處聲壓級明顯高于L-曲線法。該結果說明若在明區達到相同的期望聲壓級，GCV法中揚聲器陣列需要更高的功率，與性能仿真分析結果一致。

實驗時，使用FPGA 系統產生信號經功率放大后驅動16獨立通道揚聲器陣列，信號測量采用聲傳科技公司的CHZ-215型傳聲器頭和YG-201型前置放大器組成的寬帶測量傳聲器。測量結果如圖7(b)所示，與仿真結果(圖7(a))基本一致。在這三種方法中，未控制的聲場整體聲壓級最小，且均勻度較差，地面位置2.5 m處甚至出現較明顯的谷點。基于GCV 的正則化參數法均勻度要優于L-曲線法，該聲壓級波動范圍為8 dB。但是因為實際實驗中揚聲器功率限定，該平均聲壓級低于期望聲壓85 dB。L-曲線法的重建聲壓級整體高于GCV 法5 dB，且在期望聲壓級上下波動，波動范圍為10 dB。

為確保算法的正確性，實驗分別采集頻率為500 Hz、1500 Hz 和2000 Hz 時的聲壓數據并與1000 Hz 進行對比，如圖8 所示。可以看出相同功率下，在0～1000 Hz 范圍內L-曲線法整體聲壓級高于GCV 法5 dB 以上，并且高頻段L-曲線法均勻度以及重建誤差與GCV 法持平。另一方面L-曲線法在不同頻率下重建聲壓值皆在85～90 dB 范圍內，而GCV法波動較大。因此L-曲線法在實際應用中更具有優越性。該實驗證明，基于L-曲線法的最小二乘法聲場重建算法在功率限定的條件下，陣列能量效率較高，更符合實際需求，與性能仿真結論一致。

圖6 L-曲線法和GCV 法均勻聲場重建仿真Fig.6 Simulation of uniform sound field reconstruction using L-curve method and GCV method

圖7 三種方法的明區聲壓級分布Fig.7 Sound pressure level distribution in bright area of three methods

圖8 不同采集頻率下聲場測試結果Fig.8 Sound field test results at different acquisition frequencies

4 結論

針對聲場重建模型正則化參數選取的問題，本文將L-曲線參數選擇法引入聲場重建算法，并對比了GCV 與L-曲線法聲場重建性能。相對于GCV法，L-曲線法更具有優勢，計算L-Corner 是一個容易定義的數值問題，并且該方法很少受到相關性誤差的影響。本文以均勻度、明區能量占比以及明區重建誤差作為評價標準，對L-曲線法與GCV 法以及最小二乘法進行仿真以及實驗比較。仿真以及實驗結果顯示，L-曲線法實現了重建誤差與揚聲器驅動信號功率之間的平衡，克服了低頻段不穩定的問題，其缺點是低頻段均勻度要差于GCV法。