基于L-曲線參數(shù)優(yōu)化的均勻聲場(chǎng)重建算法

張 旭 趙洪亮 程希光 裴景春

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

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

0 引言

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

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

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

1 基于L-曲線法的均勻聲場(chǎng)重建算法

1.1 均勻聲場(chǎng)重建模型

將控制區(qū)域劃分為明區(qū)以及暗區(qū)。明區(qū)為目標(biāo)聲場(chǎng)，為達(dá)到良好的均勻聲場(chǎng)，該期望聲壓幅值設(shè)為1?；谧畲蠡刂颇繕?biāo)區(qū)域聲場(chǎng)能量的目的，設(shè)置暗區(qū)并將暗區(qū)期望聲壓設(shè)為0。聲場(chǎng)重建區(qū)域明暗區(qū)劃分如圖1所示。

圖1 聲場(chǎng)重建區(qū)域明暗區(qū)示意圖Fig.1 Schematic diagram of the light and dark area of the sound field reconstruction

在控制區(qū)域0～B3內(nèi)，明區(qū)范圍B1～B2設(shè)置M個(gè)控制點(diǎn)，暗區(qū)0～B1以及B2～B3設(shè)置N個(gè)控制點(diǎn)。明暗區(qū)通過(guò)控制點(diǎn)到陣元的距離(m=1,2,··· ,Mb)和(m=1,2,··· ,Nd)得到相應(yīng)的陣列響應(yīng)和而對(duì)于L個(gè)揚(yáng)聲器陣元權(quán)重gL，兩者關(guān)系表示如下：

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

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

本文優(yōu)化模型是在考慮最小化控制區(qū)域重建聲場(chǎng)的誤差函數(shù)的基礎(chǔ)上，約束揚(yáng)聲器權(quán)重的平方和上限，其優(yōu)化模型如下：

其中，∥·∥為2-范數(shù)，g為揚(yáng)聲器權(quán)向量，K1為其揚(yáng)聲器功率限制參數(shù)。

1.2 L-曲線選取正則化參數(shù)

為獲得揚(yáng)聲器陣列權(quán)重，需對(duì)優(yōu)化模型進(jìn)行求解。式(3)具有二次目標(biāo)和單個(gè)二次約束，因此該優(yōu)化模型為凸優(yōu)化問(wèn)題。式(3)對(duì)應(yīng)的廣義拉格朗日罰函數(shù)可以寫為

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

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

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

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

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

其中，μ′、v′和μ′′、v′′分別是μ、v對(duì)λ的一階、二階導(dǎo)數(shù)。將求得的參數(shù)值λL帶入式(5)得到揚(yáng)聲器陣列權(quán)向量gL。

1.3 聲場(chǎng)均勻控制的性能指標(biāo)

為衡量揚(yáng)聲器陣列在控制區(qū)域的均勻聲場(chǎng)重建效果，本文從以下三個(gè)參數(shù)對(duì)重建聲場(chǎng)性能進(jìn)行評(píng)價(jià)，分別為明區(qū)重建均方誤差、陣列能量效率以及聲場(chǎng)均勻度。重建性能指標(biāo)如下。

1.3.1 明區(qū)重建均方誤差

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

式(13)中，pb、Pb分別為明區(qū)重建聲場(chǎng)聲壓與明區(qū)期望聲壓；Hb為明區(qū)傳遞函數(shù)矩陣。該指標(biāo)用于衡量聲場(chǎng)重建效果的準(zhǔn)確度。

1.3.2 聲能量陣列效率

除了明區(qū)重建誤差外，仿真中還需要比較各種情況下?lián)P聲器陣列輻射到目標(biāo)區(qū)域的聲能量大小，本文使用揚(yáng)聲器陣列在明區(qū)的聲能量效率進(jìn)行評(píng)價(jià)，明區(qū)控制點(diǎn)上平均聲能量可以表示為

式(14)中，Q為目標(biāo)區(qū)域控制點(diǎn)的數(shù)量，ptm為明區(qū)中第m個(gè)控制點(diǎn)的陣列響應(yīng)，Hb為對(duì)應(yīng)的揚(yáng)聲器陣列單元到明區(qū)內(nèi)控制點(diǎn)的傳遞函數(shù)矩陣。陣列能量效率的定義為明區(qū)M個(gè)控制點(diǎn)的平均聲能量與揚(yáng)聲器陣列輸出總能量的平均聲能量的比值，其表達(dá)式如下：

該方程代表了揚(yáng)聲器陣列輻射到明區(qū)的聲能量效率，即對(duì)總能量的利用效率。方程中對(duì)控制區(qū)域控制點(diǎn)的能量取均值是為了避免控制點(diǎn)數(shù)量的變化對(duì)該參數(shù)的影響。

1.3.3 明區(qū)聲壓級(jí)均勻長(zhǎng)度和聲壓級(jí)波動(dòng)標(biāo)準(zhǔn)差

對(duì)聲場(chǎng)均勻度進(jìn)行評(píng)價(jià)時(shí)，主要采用明區(qū)的均勻長(zhǎng)度和聲壓級(jí)波動(dòng)的標(biāo)準(zhǔn)差兩個(gè)參數(shù)。明區(qū)的均勻長(zhǎng)度定義為明區(qū)聲壓幅值與均值相差少于變動(dòng)范圍3 dB 內(nèi)的距離長(zhǎng)度，該指標(biāo)可以反映出明區(qū)中均勻聲場(chǎng)所占的比例。而明區(qū)聲壓級(jí)波動(dòng)的標(biāo)準(zhǔn)差則可以反映出重建區(qū)域整體的聲場(chǎng)波動(dòng)情況。明區(qū)的均勻長(zhǎng)度rT和聲壓級(jí)波動(dòng)的標(biāo)準(zhǔn)差Sd分別表示為

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

2 數(shù)值仿真

2.1 L-曲線法的參數(shù)仿真

如圖1 所示，取揚(yáng)聲器陣元個(gè)數(shù)L=16，激勵(lì)頻率f=1000 Hz，陣元間距d=0.1 m，陣列高度H1=7 m。在控制區(qū)0～60 m 內(nèi)，設(shè)置10～50 m為明區(qū)范圍，即保持聲場(chǎng)均勻的區(qū)域，其余為暗區(qū)范圍。控制區(qū)期望聲場(chǎng)Pc=[Pb,Pd]H。Pb、Pd分別為明區(qū)期望聲壓和暗區(qū)期望聲壓?；诎祬^(qū)聲壓數(shù)值盡量小的考慮，因此設(shè)置Pd=[0,0,··· ,0N]H。明區(qū)期望聲壓表示如下：

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

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

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

圖3 λL 取值Fig.3 λLValue

2.2 重建性能指標(biāo)仿真

激勵(lì)頻率范圍設(shè)置為0～2000 Hz，為有效重現(xiàn)該范圍內(nèi)的聲場(chǎng)，選取陣元間距為d=0.1 m。根據(jù)定義的性能指標(biāo)，在0～2000 Hz 頻帶上，分別對(duì)傳統(tǒng)最小二乘法、GCV 法以及L-曲線法進(jìn)行仿真比較。

如圖4(a)、圖4(b)所示，可以明顯看出在低頻段約0～600 Hz，未正則化的最小二乘法聲場(chǎng)重建均勻度存在明顯的波動(dòng)，而高頻段聲場(chǎng)均勻度隨著頻率的變化較為連續(xù)且對(duì)應(yīng)的曲線較為平滑。陣列能量效率在0～600 Hz 頻段內(nèi)較低，波動(dòng)范圍為?160～?140 dB，而在600～1400 Hz 間，陣列效率頻率上升而提升至?20 dB?；谝陨戏治隹梢钥闯?，傳統(tǒng)的最小二乘法存在低頻段能量過(guò)低以及均勻度波動(dòng)問(wèn)題。

圖4 重建性能指標(biāo)仿真對(duì)比結(jié)果Fig.4 Reconstruction performance index simulation comparison result

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

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

3 實(shí)驗(yàn)比較

本節(jié)將揚(yáng)聲器陣列未控制的聲場(chǎng)效果與GCV法及L-曲線法的重建聲場(chǎng)進(jìn)行比較。本實(shí)驗(yàn)將明區(qū)聲場(chǎng)縮小為1～5 m，陣元中心距離地面為2 m，揚(yáng)聲器陣列數(shù)量為16個(gè)，陣元間距為0.6 m，測(cè)量間距為0.2 m，如圖5所示。期望聲壓設(shè)為85 dB，并將三種方法的揚(yáng)聲器功率均設(shè)置為以L-曲線功率為標(biāo)準(zhǔn)進(jìn)行計(jì)算。實(shí)驗(yàn)以及仿真頻率均取1000 Hz，在此頻率下，L-曲線法取正則化參數(shù)值為λL=0.4480，GCV法取參數(shù)值為λG=2.66×10?5。

圖5 實(shí)驗(yàn)裝置圖Fig.5 Experimental setup diagram

首先對(duì)L-曲線法及GCV 法進(jìn)行均勻聲場(chǎng)重建仿真。如圖6所示，該聲場(chǎng)左側(cè)聲壓最大值處為揚(yáng)聲器陣列，陣列根據(jù)求得的揚(yáng)聲器權(quán)重偏轉(zhuǎn)一定角度，在明區(qū)范圍內(nèi)進(jìn)行均勻聲場(chǎng)重建。由圖6 的聲壓級(jí)分布圖可以看出，L-曲線法在明區(qū)重建區(qū)域誤差以及均勻度略差于GCV 法。然而，GCV 法的陣列處聲壓級(jí)明顯高于L-曲線法。該結(jié)果說(shuō)明若在明區(qū)達(dá)到相同的期望聲壓級(jí)，GCV法中揚(yáng)聲器陣列需要更高的功率，與性能仿真分析結(jié)果一致。

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

為確保算法的正確性，實(shí)驗(yàn)分別采集頻率為500 Hz、1500 Hz 和2000 Hz 時(shí)的聲壓數(shù)據(jù)并與1000 Hz 進(jìn)行對(duì)比，如圖8 所示?？梢钥闯鱿嗤β氏拢?～1000 Hz 范圍內(nèi)L-曲線法整體聲壓級(jí)高于GCV 法5 dB 以上，并且高頻段L-曲線法均勻度以及重建誤差與GCV 法持平。另一方面L-曲線法在不同頻率下重建聲壓值皆在85～90 dB 范圍內(nèi)，而GCV法波動(dòng)較大。因此L-曲線法在實(shí)際應(yīng)用中更具有優(yōu)越性。該實(shí)驗(yàn)證明，基于L-曲線法的最小二乘法聲場(chǎng)重建算法在功率限定的條件下，陣列能量效率較高，更符合實(shí)際需求，與性能仿真結(jié)論一致。

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

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

圖8 不同采集頻率下聲場(chǎng)測(cè)試結(jié)果Fig.8 Sound field test results at different acquisition frequencies

4 結(jié)論

針對(duì)聲場(chǎng)重建模型正則化參數(shù)選取的問(wèn)題，本文將L-曲線參數(shù)選擇法引入聲場(chǎng)重建算法，并對(duì)比了GCV 與L-曲線法聲場(chǎng)重建性能。相對(duì)于GCV法，L-曲線法更具有優(yōu)勢(shì)，計(jì)算L-Corner 是一個(gè)容易定義的數(shù)值問(wèn)題，并且該方法很少受到相關(guān)性誤差的影響。本文以均勻度、明區(qū)能量占比以及明區(qū)重建誤差作為評(píng)價(jià)標(biāo)準(zhǔn)，對(duì)L-曲線法與GCV 法以及最小二乘法進(jìn)行仿真以及實(shí)驗(yàn)比較。仿真以及實(shí)驗(yàn)結(jié)果顯示，L-曲線法實(shí)現(xiàn)了重建誤差與揚(yáng)聲器驅(qū)動(dòng)信號(hào)功率之間的平衡，克服了低頻段不穩(wěn)定的問(wèn)題，其缺點(diǎn)是低頻段均勻度要差于GCV法。