霍茹靜，魏明洋，許志勇

(1.南京理工大學電子工程與光電技術(shù)學院，江蘇南京 210094；2.中國科學院聲學研究所，北京 100190)

0 引 言

隨著國內(nèi)民用車輛數(shù)量的增加，車輛私自改裝、大型卡車超載、車輛老化不報廢等問題日益突出，其帶來的噪聲污染、妨礙安全行車等問題亟待解決。這些問題車輛在道路上行駛時都會產(chǎn)生異常聲音，所以現(xiàn)階段的解決方案是通過異常聲檢測[1]與方位估計(Direction of Arrival, DOA)算法[2]確定問題車輛的位置和違法情況。由于異常聲源往往處于移動狀態(tài)，對位置信息的確定要求具有實時性；同時通過對這類異常聲進行分析統(tǒng)計發(fā)現(xiàn)，異常聲的聲源級高、傳播距離遠、聲信號的能量集中于低頻段(500～1 000 Hz)且?guī)捿^寬[3]，但是低頻信號的高精度DOA估計需要使用大尺寸傳聲器陣列對聲音進行采集。在實際應用中，需要結(jié)合上述特點選擇合適的DOA算法與傳聲器陣型，確保DOA估計的高精度和低運算復雜度。

20世紀70年代末，多重信號分類(Multiple Signal Classification, MUSIC)[4]算法的出現(xiàn)開啟了超分辨DOA估計算法的篇章。該算法通過譜峰搜索得到目標方位。當對精度要求較高時，譜峰搜索步長很小，這會導致算法的運算量很大。旋轉(zhuǎn)不變子空間技術(shù)(Estimation of Signal Parameters via Rotational Invariance Techniques, ESPRIT)[5-6]算法相比MUSIC算法無需進行譜峰搜索，計算速度得到較大提高且DOA估計誤差小，因而得到了廣泛應用。

分布式陣列[7-8]是由空間散布的多個子陣列構(gòu)成的陣列系統(tǒng)，其靈活的空間布置特性使得分布式陣列廣泛應用于復雜安裝環(huán)境的場景，子陣的分開布放使得分布式陣列具有很大的陣列尺寸，從而具有很好的低頻信號DOA估計性能。分布式陣列的陣元總數(shù)多，陣列具有很強的魯棒性以及很高的陣增益。這些優(yōu)勢使得基于分布式陣列的DOA估計算法研究具有重要意義。但是大尺寸陣列的陣元間距大于信號半波長，不滿足空間采樣定理，在對信號進行DOA估計時會同時出現(xiàn)多個估計值，即測角模糊。為了得到無模糊且高精度的DOA估計結(jié)果，需要針對陣型采取解模糊[9-11]策略。

本文提出了一種用于公路異常聲監(jiān)測系統(tǒng)的分布式陣列DOA估計算法，利用ESPRIT算法無需進行譜峰搜索的特點降低了計算復雜度。選取三個均勻矩形陣(Uniform Rectangular Array, URA)呈三角形布放在公路監(jiān)控桿上，陣列整體孔徑大且具有空間平移不變性。該算法利用寬帶聚焦[12-13]方式將經(jīng)典ESPRIT算法的適用范圍從窄帶信號拓展至寬帶信號，并利用陣列兩種尺度下的方位信息旋轉(zhuǎn)不變特性估計聲源方位，再根據(jù)最小誤差準則解模糊法實現(xiàn)公路異常聲源的高精度DOA估計。本文分析了分布式陣列較單個陣列對低頻信號DOA估計精度的提升，研究了子陣間基線距離對DOA估計精度的影響。仿真結(jié)果驗證了本文算法的有效性。

1 分布式陣列數(shù)學模型

如圖1所示，子陣呈三角形分布的傳聲器陣列擺放在監(jiān)控桿上，假設(shè)陣列分布在yOz平面上，坐標原點O位于子陣S1的第一個傳聲器處，子陣S1，S2，S3為孔徑相等的URA，均包含M2個傳聲器，單個子陣的傳聲器間距d≤λ/2，λ為聲源信號分析頻段中最大頻率對應的波長。在y軸方向上兩個子陣的基線距離為Dy1，且滿足Dy1?Md，在z軸方向上的子陣S3相對于子陣S1的基線距離在y軸和z軸方向上的投影分別為Dy2和Dz，且滿足Dy2?Md,Dz?Md。

圖1 位于交通監(jiān)控桿上的分布式傳聲器陣列示意圖Fig.1 Schematic diagram of the distributed microphone array on traffic monitoring pole

假設(shè)一個平面波寬帶聲源入射到該分布式陣列上，記來波方向在xOy平面上的投影與x軸正半軸的夾角θ為方位角，來波方向與z軸正半軸的夾角φ為俯仰角，y軸與z軸上的方向余弦分別為u=sinθsinφ，v=cosφ，方位角與俯仰角的取值范圍分別為θ∈[0，π2]，φ∈[0，π2)。

由于陣列接收的公路聲源信號通常是非平穩(wěn)的，因此需要先對接收信號進行分幀，然后逐幀變換到頻率域。定義信號能量主要分布頻段的帶寬為B，將該頻段劃分為J個頻點，可以得到第i幀的寬帶信號模型為

式中：Xi(fj)、Si(fj)、Ni(fj)分別為陣列接收數(shù)據(jù)、聲源信號、噪聲在fj頻點處的頻譜向量，j=1，2，…，J。Ai(fj)為陣列流形矩陣，本文假定聲源個數(shù)為1，則矩陣退化為列向量a(u，v，fj)。令子陣S1的第一行與第一列的導向矢量分別為

則子陣S1，S2，S3的導向矢量分別為

可得：

式中：?表示克羅內(nèi)克(Kronecker)積。

在第i幀的頻點fj下，X(fj)的協(xié)方差矩陣為

其中，信號能量σS2(fj)=Si(fj)Si*(fj)，Rn，i(fj)是噪聲協(xié)方差矩陣。假設(shè)聲源信號在連續(xù)2K+1幀時間內(nèi)是廣義平穩(wěn)的，則第i幀協(xié)方差矩陣可用前后K幀信號協(xié)方差矩陣的統(tǒng)計平均來估計，即：

2 分布式二維DOA估計算法

2.1 寬帶聚焦算法

公路車輛異常聲信號為寬帶信號，在利用ESPRIT算法進行DOA估計時，需要先將頻帶內(nèi)各頻點的信號子空間聚焦到參考頻點，得到基于參考頻點的樣本協(xié)方差矩陣，進而根據(jù)子空間估計波達方向。本文采用雙邊相關(guān)變換(Two-sided Correlation Transform, TCT)聚焦算法[14]，該算法在理想條件下滿足：

而實際情況中，樣本協(xié)方差矩陣中包含噪聲，聚焦變換矩陣要滿足歸一化約束，即要滿足：

其中：T(fj)=U(f0)UH(fj)，U(f0)和U(fj)分別是和的最大奇異值對應的奇異向量。

2.2 二維經(jīng)典ESPRIT算法

應用ESPRIT算法的核心思想是構(gòu)造兩個相同的子陣，一子陣的導向向量乘以關(guān)于方位信息的旋轉(zhuǎn)因子可以得到另一子陣的導向向量(即空間平移不變)，二維ESPRIT算法要求在兩個方向上構(gòu)造兩組具有空間平移不變性的子陣。

分布式均勻矩形陣布放示意圖如圖2所示。在圖2中，在y軸方向選取子陣S1左邊M-1列傳聲器構(gòu)成子陣SL，右邊M-1列傳聲器構(gòu)成子陣SR，SL與SR之間的位移量為d。

圖2 分布式均勻矩形陣布放示意圖Fig.2 Layout diagram of the distributed URA

則陣列SL與SR的導向矢量滿足關(guān)系

其中，JSL與JSR分別表示陣列SL與SR的選擇矩陣，JSL=I2×3?(IM?[IM-10(M-1)×1])，I為單位矩陣，JSR=I2×3?(IM?[0(M-1)×1IM-1])，Φ=ej2πduλ表示陣列SL與SR的平移關(guān)系，其中波長λ=c/f0，c為聲速，f0為聚焦頻點。

聚焦后的信號協(xié)方差矩陣R(f0)進行特征值分解：

其中，ΣS和ΣN分別表示R(f0)的最大特征值和小特征值組成的對角陣，最大特征值對應的特征向量為信號子空間US。由于導向矢量a(u，v)與子空間JSLUS以及子空間JSRUS張成的空間相同，即span{a(u，v)}=span{JSLUS}=span{JSRUS}，span{}·表示生成向量張成的空間，則信號子空間的旋轉(zhuǎn)不變特性可以表示為

根據(jù)最小二乘法得到旋轉(zhuǎn)矩陣：

由于信號源個數(shù)為1，Ψ矩陣退化為標量Ψ。定義空間角頻率β=angle(Ψ)，angle(·)表示求相位角。則y軸的方向余弦估計結(jié)果為

同理在z軸方向，選取子陣S1和S3下方M-1行傳聲器構(gòu)成子陣SD，上方M-1行構(gòu)成子陣SU，SD與SU之間的位移量為d，可以求得z軸的方向余弦估計結(jié)果v。

2.3 二維雙尺度ESPRIT算法

采用圖2中的分布式傳聲器陣列，利用子陣內(nèi)相鄰陣元的間距作為陣列平移不變尺度可以得到精度較低但無模糊的方向余弦估計值，稱為粗估計；利用兩兩子陣間的間距作為陣列平移不變尺度可以得到精度高但有模糊的估計值，稱為精估計。雙尺度ESPRIT算法先對陣列接收信號進行粗估計和精估計，再依據(jù)最小誤差準則并結(jié)合粗估計結(jié)果對精估計結(jié)果解模糊，得到最終DOA估計結(jié)果。

2.3.1 粗估計與精估計

根據(jù)定義，由2.2節(jié)所得的方向余弦估計值u和v可記為方向余弦粗估計uc和vc。類似地，精估計的空間平移不變性可表示為

信號子空間的旋轉(zhuǎn)不變特性可表示為

其中，JSf1與JSf2表示精估計的選擇矩陣，Ψf為精估計的平移關(guān)系。由于信號源個數(shù)為1，Ψf退化為標量Ψf，精估計相位為βf=angle(Ψf)。

定義子陣S1，S2，S3的選擇矩陣分別為JS1= [IMM,02MM×MM]，JS2= [0MM,IMM,0MM]，JS3= [0MM,0MM,IMM]，0表示元素全為0的矩陣，則利用兩兩子陣進行精估計的選擇矩陣為(JSf1，JSf2)=(JS1，JS2)，(JSf1，JSf2)=(JS1，JS3)，(JSf1，JSf2)=(JS2，JS3)，對 應 的Φf分 別 為Φ12=ej2πDy1uλ，Φ13=ej2π(Dy2u+Dzv)λ，Φ23=ej2π((Dy1+Dy2)u+Dzv)λ，將精估計選擇矩陣依次代入式(21)求得角頻率βf12、βf13、βf23。

2.3.2 最小誤差準則解模糊

由于Dy1、Dy2、Dz均遠大于d，在估計相位βf時會產(chǎn)生以2π為周期的測角模糊，則兩兩子陣間角頻率βf12、βf13、βf23滿足關(guān)系：

其中：n1，n2，n3∈Z，uf與vf為存在模糊的方向余弦精估計。聯(lián)立式(22)中前兩式即可求解無模糊的方向余弦精估計udf與vdf。

由圖3所示在平面(uf，vf)上，式(22(a))表示平行于vf軸的均勻直線簇，(22(b))表示斜率相等而截距不同的均勻直線簇，方向余弦精估計為單位圓內(nèi)兩直線簇的交點，粗估計結(jié)果在該平面上表示參考點(uc，vc)。對方向余弦精估計解模糊，即根據(jù)最小誤差準則，找到距參考點最近的直線交點。為實現(xiàn)正確解模糊，參考點應滿足：

圖3 解模糊法示意圖Fig.3 Schematic diagram of defuzzification method

若直接按式(23)的約束遍歷所有直線交點，計算復雜度較高。為此，本文提出一種結(jié)合幾何意義的解模糊方法，首先確定周期模糊數(shù)n1，即找到距參考點最近的斜率為0的直線L1：

此時，參考點(uc，vc)在直線L1上的映射為(udf，vc)。

然后，確定周期模糊數(shù)n2，即找到距映射點最近的斜率不為0的直線L2：

最后得到聲源DOA估計的方位角φDOA和俯仰角θDOA：

2.4 算法流程

綜上所述，基于雙尺度ESPRIT的寬帶聲源DOA估計算法流程總結(jié)如下：

(1) 將信號頻帶內(nèi)各頻點處的信號子空間聚焦到參考頻點f0，并對聚焦后的協(xié)方差矩陣R(f0)進行特征值分解，得到信號子空間。

(2) 根據(jù)式(16)分別構(gòu)造沿y軸和z軸粗估計的旋轉(zhuǎn)不變性等式，求出方向余弦粗估計uc和vc。

(3) 對于S1，S2，S3中的兩兩子陣，根據(jù)式(20)構(gòu)造精估計的旋轉(zhuǎn)不變性等式，得到存在周期模糊的相位βf12、βf13、βf23。

(4) 由式(22)結(jié)合幾何意義，根據(jù)最小誤差準則求得無模糊的方向余弦精估計udf和vdf。

(5) 由式(28)求出異常聲源DOA估計的方位角和俯仰角φDOA和θDOA。

3 仿真分析

設(shè)M=4，三個子陣內(nèi)傳聲器間距d=0.01 m，分布式陣列傳聲器總數(shù)為3M2=48。子陣間基線距離滿足0.04 m＜Dy1+Dy2＜4 m，0.04 m＜Dz＜3.5 m。模擬車輛異常聲寬帶信號源的入射方向(θ，φ)=(30°，60°)，信噪比RSN=10 dB。采用聯(lián)合角度估計與方向余弦的均方根誤差(Root Mean Square Error,RMSE)作為衡量算法的估計精度性能，定義為

其中，K表示蒙特卡洛(Monte Carlo)實驗次數(shù)本文中K=5 000，與為第k次實驗信號的估計值。

3.1 分布式DOA估計的性能比較與驗證

本次實驗設(shè)置分布式陣列的基線距離分別為Dy2=0.3 m，Dy1=1.5 m，Dz=1.5 m，此外，設(shè)置陣元數(shù)和陣元間距與分布式陣列相同的6行8列的URA作為對比陣列。圖4給出了兩個陣列基于ESPRIT算法的聯(lián)合DOA估計誤差在不同頻率范圍(1/3倍頻程帶寬)的變化。圖4中的RMSE根據(jù)式(29)計算，橫坐標為1/3倍頻程帶寬對應的中心頻率。

圖4 各頻帶下分布式陣列的聯(lián)合DOA估計均方根誤差Fig.4 RMSE of joint DOA estimation of distributed array in each frequency band

由圖4可以看出，DOA估計誤差隨信號頻率的增大而減小，分布式陣列采用雙尺度ESPRIT算法的估計精度比URA提高了約10倍，而僅采用粗估計尺度經(jīng)典ESPRIT算法的分布式陣列估計精度比URA差，這是由于粗估計只利用了(24，24)個陣元估計u和v，而URA使用了(42，40)個陣元。該實驗驗證了本文提出算法以及通過增大陣列孔徑的方式提高DOA估計精度的可行性。

3.2 估計誤差隨基線距離變化情況

下面探究分布式陣列ESPRIT算法的估計精度與基線距離的關(guān)系。仿真使用的低頻寬帶聲源信號的帶寬為562～708 Hz。

(1) 仿真假設(shè)Dy2=0.3 m，圖5給出了方向余弦u的估計誤差隨Dy1/d的變化情況。圖5中的RMSE根據(jù)式(30)計算。由圖5可以看出，隨著Dy1/d的增大，在一定范圍內(nèi)，雙尺度ESPRIT算法的方向余弦u的估計精度逐漸提升，但當d＞1.9 m、Dy1/d＞190時，估計精度變差，此時的基線距離稱為基線模糊門限。經(jīng)典ESPRIT算法只利用了相鄰陣元間距，基線距離Dy1對方向余弦u的估計精度影響不大。

圖5 方向余弦u估計精度與Dy1/d的關(guān)系Fig.5 Relationship between accuracy of direction cosine estimation u and baseline distance Dy1

(2) 令Dy1=1.9 m，圖6給出了在不同基線距離Dy2下，方向余弦v的估計精度隨基線距離Dz的變化情況。圖6中的RMSE根據(jù)式(31)計算。由圖6可以看出，在基線模糊門限范圍內(nèi)，Dy2=0時的估計精度高于Dy2=0.3 m時的估計精度。

圖6 方向余弦v估計精度與基線距離Dz/d的關(guān)系Fig.6 Relationship between accuracy of direction cosine estimation v and baseline distance Dz

由于解模糊以粗估計值為參考，當基線距離Dy1與Dz增大到一定程度，圖3中的直線簇變密集，當粗估計值無法滿足式(23)中的條件時，會導致解模糊錯誤率升高，DOA估計誤差增大。當Dy2=0，圖3中的直線簇相互垂直，式(27)中v不受u的誤差影響，此時的DOA估計精度最高。在實際應用中，可根據(jù)公路監(jiān)控桿的條件設(shè)計分布式陣型。

4 結(jié) 論

為解決公路異常聲源的高精度DOA估計問題，并滿足低計算復雜度、布陣靈活的需求，本文提出了一種適用于均勻矩形子陣呈三角形分布的分布式傳聲器陣列二維DOA估計算法。本文對信號采用寬頻段聚焦和ESPRIT算法估計DOA，避免了逐頻點估計和二維譜峰搜索造成的復雜計算，算法的雙尺度策略與最小誤差準則解模糊法能利用擴展孔徑的優(yōu)勢，提高低頻聲源的DOA估計精度，且適用于任意三個相同均勻矩形子陣組成的分布式陣列，具有較好的工程應用前景。