單穎春，劉耀光，劉獻(xiàn)棟，何 田

(北京航空航天大學(xué) 交通科學(xué)與工程學(xué)院，北京 100191)

采用陣列信號(hào)處理方法進(jìn)行聲源識(shí)別時(shí)，如果入射聲信號(hào)相干，許多傳統(tǒng)的自適應(yīng)陣列將出現(xiàn)聲源信號(hào)相互抵消的現(xiàn)象［1］，難以實(shí)現(xiàn)對(duì)聲源位置的正確識(shí)別。為了解決這一問題，許多學(xué)者提出了其他的陣列信號(hào)處理算法，其中應(yīng)用最為廣泛的是由 Evans［2］、Shan［3］提出的的空間平滑法，該方法能重建信號(hào)協(xié)方差矩陣的秩。Chen［4］將一維空間平滑法發(fā)展到二維。但是，以上的空間平滑法只適用于陣列是規(guī)則矩陣，信號(hào)來自于遠(yuǎn)場的情況。為此，F(xiàn)riedlander［5］結(jié)合矩陣插值技術(shù)使空間平滑法能用于任意形狀陣列，在此基礎(chǔ)上Yang［6］對(duì)一維近場相干信號(hào)的波達(dá)方向進(jìn)行估計(jì)，實(shí)現(xiàn)了一維近場聲信號(hào)的識(shí)別。

本文應(yīng)用矩形陣列采集聲信號(hào)，對(duì)二維近場相干聲源識(shí)別進(jìn)行了研究首先通過插值矩陣技術(shù)把實(shí)際陣列接收到的近場信號(hào)轉(zhuǎn)換為形式上的遠(yuǎn)場信號(hào)，然后在遠(yuǎn)場信號(hào)的基礎(chǔ)上通過二維空間平滑法求出陣列的平滑協(xié)方差矩陣，再針對(duì)平滑協(xié)方差矩陣，使陣列信號(hào)輸出功率最小，從而獲得給定方向向量對(duì)應(yīng)的聲信號(hào)功率，實(shí)現(xiàn)聲源的識(shí)別。

1 聲源信號(hào)模型

假設(shè)聲源平面上有K個(gè)點(diǎn)聲源，測試平面上用均勻分布的矩形陣列對(duì)信號(hào)進(jìn)行采集。陣列為M×N維，各傳聲器相同且忽略傳聲器的方向性帶來的影響，x，y兩方向上傳聲器的間隔分別為dx，dy。聲源與傳聲器陣列的相對(duì)位置見圖1。

根據(jù)Mailloux［7］的經(jīng)驗(yàn)公式，當(dāng)點(diǎn)聲源與陣列原點(diǎn)之間的距離r滿足r≤2L2/λ時(shí)聲源屬于近場聲源，上式中L為陣列的最大尺寸，λ為聲源的波長。此時(shí)聲源的傳播規(guī)律應(yīng)按球面波傳播進(jìn)行分析。如r≥2L2/λ則聲源屬于遠(yuǎn)場，傳播規(guī)律可按平面波傳播規(guī)律進(jìn)行分析。

當(dāng)K個(gè)聲信號(hào)從近場入射時(shí)，則:

第(m，n)個(gè)傳聲器接收到的信號(hào)xmn為:

式中:rk(m，n)是第k個(gè)聲源到第(m，n)個(gè)傳聲器之間的距離;sk(t)是第k個(gè)入射點(diǎn)聲源隨時(shí)間的變化規(guī)律;f是聲源的頻率;c為聲波在介質(zhì)中傳播的速度;nmin(t)是各通道中方差為σ2的白噪聲，各通道噪聲互不相關(guān)，由傳聲器自身和環(huán)境噪聲引起。

把各傳聲器接收到的信號(hào)表示成向量形式，見公式(2):

式中，

其中，Xn(t)是各傳聲器接受到的信號(hào)，S(t)為聲源信號(hào)向量;An為近場方向矩陣;ak是第k個(gè)聲源所對(duì)應(yīng)的方向向量。N(t)是各通道噪聲向量。

當(dāng)K個(gè)聲信號(hào)從遠(yuǎn)場入射時(shí)，傳聲器接收到的信號(hào)表示成向量形式為:

式中:S(t)、N(t)的含義、形式與公式(2)相同。Xf(t)表示聲源為遠(yuǎn)場時(shí)傳聲器接受到的信號(hào)，其形式與公式(2)中的Xn(t)相同。

其中，θk，φk分別為第k個(gè)聲源的方位角和仰角。

2 近場相干聲源的定位方法

空間平滑法只能對(duì)陣列接受到的遠(yuǎn)場聲信號(hào)進(jìn)行數(shù)據(jù)處理以實(shí)現(xiàn)聲源的識(shí)別［6］，對(duì)于近場下測得的數(shù)據(jù)必須通過矩陣插值技術(shù)，使近場方向矩陣An轉(zhuǎn)化為遠(yuǎn)場方向矩陣Af。在此基礎(chǔ)上再通過二維空間平滑法把相干信號(hào)的信號(hào)協(xié)方差矩陣中非對(duì)角線項(xiàng)元素置零，然后估計(jì)信號(hào)在給定方向下的功率。

2.1 插值矩陣技術(shù)

插值矩陣技術(shù)［5］的基本思想是將一個(gè)矩陣的每一個(gè)列向量線性表示另一個(gè)矩陣的列向量，一般通過最小二乘法找出兩個(gè)矩陣之間的插值矩陣(轉(zhuǎn)換矩陣)。假定平面的K個(gè)相干點(diǎn)聲源集中分布在某一區(qū)域內(nèi)，這一區(qū)域如圖2所示，聲源相對(duì)于測試陣列的方位角、仰角的范圍分別為:θrange=［θmin，θmax］，φrange=［φmin，φmax］。例如，要確定一輛汽車側(cè)面聲源的位置，那么方向角、仰角的范圍將由傳聲器陣列的大小、安裝位置，傳聲器陣列與汽車之間距離，汽車側(cè)面的大小來決定。

為了得到近場信號(hào)與遠(yuǎn)場信號(hào)之間的轉(zhuǎn)換矩陣，分別以間隔Δθ，Δφ對(duì)方位角和仰角進(jìn)行劃分，如圖3所示。

現(xiàn)假定以上區(qū)域每一交點(diǎn)上都存在一點(diǎn)聲源，則傳聲器陣列接收到信號(hào)為:

式(5)中:

其中:An-virtual表示假定每一交點(diǎn)上存在一聲源時(shí)的近場方向矩陣，Sk-virtual表示第k個(gè)交點(diǎn)上假定的聲源，Q為聲源區(qū)域上交點(diǎn)的總數(shù)。

對(duì)于同樣的信號(hào)Svirtual(t)，如果其來自于遠(yuǎn)場，則傳聲器陣列接收到的信號(hào)為:

其中，Af-virtual為假定遠(yuǎn)場聲源的方向矩陣，Af-virtual=［b1-virtual，b2-virtual，…，bQ-virtual］。

假設(shè)存在常矩陣T，使THAn-virtual=Af-virtual，則通過常矩陣T可以使實(shí)際信號(hào)與陣列間的近場方向矩陣An轉(zhuǎn)化為遠(yuǎn)場方向矩陣Af。實(shí)際上常矩陣T不一定存在，可以根據(jù)minTHAn-virtual-Af-virtualF給出常矩陣T的最小二乘解［6］:

其中:F表示矩陣的Frobenius范數(shù)。最小二乘法的性能與An-virtual、Af-virtual兩矩陣直接相關(guān)，矩陣An-virtual和Af-virtual是由角度間隔 Δθ，Δφ、重建區(qū)域大小、傳聲器陣列決定的。Δθ，Δφ由矩陣An-virtual到Af-virtual的精度決定，在 Δθ=Δφ =0.1°時(shí)可以得到滿意的結(jié)果。另外，在傳聲器陣列參數(shù)(位置、個(gè)數(shù)、形狀)已知的情況下，如果重建區(qū)域過大，將導(dǎo)致矩陣THAn-virtual、Af-virtual兩者誤差的平方和較大，這說明THAn-virtual不能很好地近似Af-virtual，可以通過減少重建區(qū)域大小來解決此問題［5］。

對(duì)式(2)左乘TH，把實(shí)際在近場中接收到的信號(hào)Xn(t)轉(zhuǎn)換為遠(yuǎn)場信號(hào)Xf(t)，即:

以上步驟為進(jìn)一步應(yīng)用空間平滑法去除信號(hào)的相關(guān)提供了必要的條件。

2.2 二維空間平滑法

二維空間平滑法首先把規(guī)則M×N維矩形陣分割成(M-Ms+1)(N-Ns+1)個(gè)互相重疊的子矩陣，每個(gè)子矩陣的大小為Ms×Ns，如圖4所示，然后分別求出每個(gè)子矩陣的協(xié)方差矩陣后相加求平均，得到二維空間平滑協(xié)方差矩陣，為了討論的簡單，本文首先忽略各通道中原有的噪聲分量。

圖4 二維空間平滑法的矩陣分割Fig.4 Division of array on 2-D spatial smoothing

第(m，n)個(gè)子陣列接收到的信號(hào)為:

其中，As-f為子陣列(1，1)的方向矩陣:

第(m，n)個(gè)子矩陣的協(xié)方差矩陣為:

其中，RSS=E［S(t)S(t)H］表示信號(hào)的協(xié)方差矩陣，當(dāng)信號(hào)互不相關(guān)時(shí)是一個(gè)對(duì)角矩陣，當(dāng)信號(hào)相干時(shí)是非對(duì)角矩陣。

定義陣列的二維空間平滑協(xié)方差矩陣為:

其中:

是信號(hào)的平滑協(xié)方差矩陣。

當(dāng)子陣列數(shù)量足夠多時(shí)，可以得到以下近似形式:

式中，RSS(i，i)=E［si(t)si(t)H］=為聲源si(t)的功率。

2.3 聲源信號(hào)的功率預(yù)測

空間平滑法后得到空間平滑協(xié)方差矩陣，在此基礎(chǔ)上，根據(jù)總功率P=WHW最小以及對(duì)應(yīng)方向向量受約束WHb=1的原則即可求出對(duì)應(yīng)方向向量的信號(hào)功率。

［8］，對(duì)應(yīng)方向向量信號(hào)的功率為:

實(shí)際應(yīng)用中R-1不存在，可以對(duì)(R+I)求逆，I為適當(dāng)?shù)脑肼晠f(xié)方差矩陣。由于實(shí)際的近場矩陣已經(jīng)轉(zhuǎn)換為遠(yuǎn)場矩陣，因而每一個(gè)遠(yuǎn)場方向向量b對(duì)應(yīng)一個(gè)近場方向向量(即對(duì)應(yīng)聲源面上的一點(diǎn))，遠(yuǎn)場方向向量b在一定方位角與仰角范圍內(nèi)掃描后就能覆蓋對(duì)應(yīng)聲源平面上一定區(qū)域，多次重復(fù)使用式(14)即得到這一區(qū)域上對(duì)應(yīng)功率分布，進(jìn)而得出聲源的分布。

3 仿真分析

為了驗(yàn)證該方法的正確性，進(jìn)行以下的仿真:三個(gè)入射點(diǎn)聲源位于Z0=1的平面上，均為1 000 Hz，且相互之間有固定相差，三個(gè)點(diǎn)聲源均位于離陣列1 m的XOY平面上，入射方位角、仰角分別為(35°，70°)，(20°，50°)，(60°，50°);對(duì)應(yīng)X、Y坐標(biāo)分別為(0.3，0.21)，(0.65，0.29)，(0.42，0.62)，陣列為 9 ×9 并且是均布的矩形傳聲器陣列，傳聲器間隔在X，Y方向上都為0.15米。

首先，采用傳統(tǒng)的延時(shí)累加波束形成，識(shí)別結(jié)果如圖5 所示，圖5 中(0.65，0.29)，(0.42，0.62)兩處的聲源沒能識(shí)別出來。其次，采用本文提出算法，聲源區(qū)域范圍取 θrange=［20°，60°］，φrange=［50°，70°］，Δθ= Δφ =0.1°，識(shí)別結(jié)果如圖6所示，準(zhǔn)確地識(shí)別出聲源位置。

由以上兩組算例可以看出，本文提出的算法能在聲源相干的環(huán)境下正常工作，且在聲源識(shí)別的精度上比傳統(tǒng)的波束形成有較大的提高。影響基于插值矩陣技術(shù)與二維空間平滑法相干聲源識(shí)別方法分辨率的因素包括聲源頻率、陣列大小等。為了比較該方法與傳統(tǒng)波束形成分辨率隨聲源頻率的變化趨勢，進(jìn)行以下仿真:聲源頻率由200 Hz取到4 000 Hz，分別計(jì)算它們的主瓣寬度，這里主瓣寬度表示只存在一個(gè)聲源時(shí)識(shí)別結(jié)果在-3 dB處曲線的最大直徑，聲源位于(0.442 3，0.371 1)，計(jì)算結(jié)果如圖7 所示，點(diǎn)劃線是傳統(tǒng)波束形成主瓣寬度隨聲源頻率的變化曲線，實(shí)線是該方法主瓣寬度隨聲源頻率的變化曲線。由圖7可以看出，盡管在某些頻率上該方法的主瓣寬度比傳統(tǒng)波束形成要大，但總的來說，該方法的主瓣寬度還是要比傳統(tǒng)波束形成小很多。

圖5 傳統(tǒng)波束形成識(shí)別結(jié)果Fig.5 Result based on traditional beamforming

圖6 基于插值矩陣技術(shù)與二維空間平滑法識(shí)別結(jié)果Fig.6 Result based on interpolated array technique and two-dimensional spatial smoothing

圖7 基于插值矩陣技術(shù)與二維空間平滑法聲源識(shí)別與傳統(tǒng)波束形成主瓣寬度隨聲源頻率的變化曲線Fig.7 The curve of main lobe width with frequency based on two methods

最后，基于插值矩陣技術(shù)與二維空間平滑法聲源識(shí)別方法要對(duì)方向角θ和仰角φ兩個(gè)變量進(jìn)行搜索，搜索八萬個(gè)方向的時(shí)間大約四分鐘。近場相干聲源通過插值矩陣技術(shù)轉(zhuǎn)換為遠(yuǎn)場相干聲源再經(jīng)過空間平滑法去相關(guān)后，如果結(jié)合超分辨率波達(dá)方向估計(jì)會(huì)使識(shí)別精度進(jìn)一步提高。

4 結(jié)論

本文首先通過插值矩陣技術(shù)把實(shí)際中接收到的近場數(shù)據(jù)轉(zhuǎn)換為形式上的遠(yuǎn)場數(shù)據(jù)，遠(yuǎn)場數(shù)據(jù)是使用二維空間平滑法的必要條件，然后在此基礎(chǔ)上通過二維空間平滑法求出陣列的平滑協(xié)方差矩陣，去除了信號(hào)之間的相關(guān)性，再由平滑協(xié)方差矩陣求出對(duì)應(yīng)某個(gè)方向向量信號(hào)的功率。通過算例證明，以上方法能對(duì)二維近場相干聲源進(jìn)行有效識(shí)別，與傳統(tǒng)波束形成法相比可以獲得更高的識(shí)別精度。在研究過程中，忽略了通道中原有的由傳聲器自身和環(huán)境引起的噪聲。

參考文獻(xiàn)

［1］ Widrow B，Duvall K M，Gooch R P，et al.Signal cancellation phenomena in adaptive antennas:causes and cures［J］.IEEE transactions on antennas and propagation，1982，30(3):469-478.

［2］Evans J E，Johnson J R，Sun D F.Application of advanced signal processing techniques to angle of arrival estimation in ATC navigation and surveillance system［R］.Technical report582，MIT:Lincoln Laboratory，1982.

［3］ Shan T J，Wax M，Kailath T.On spatial smoothing for estimation of coherent signals［J］.IEEE trans.Acoustics，Speech and Signal processing，1985，33(4):806-811.

［4］ Chen Y M.On spatial smoothing for two-dimensional directionof-arrival estimation of cohernt signals［J］.IEEE transaction on signal processing，1997，45(7):1689-1696.

［5］ Friedlander B，Weiss A J.Direction finding using spatial smoothing with interpolated arrays［J］.IEEE Transactions on Aerospace and electronic systems，1992，28(2):574-587.

［6］ Yang D S，Shi J，Liu B S.Doa estimation for the near-field correlated sources with interpolated array technique［C］.IEEE Conference on industrial electronics and applications，2009，3011-3015.

［7］ Mailloux R J.Phased array antenna handbook［M］.Norwood:Artech House Publishers，2005.

［8］Shan T J，Kailath T.Adaptive beamforming for coherent signals and interference［J］.IEEE Trans.Acoustics，Speech and Signal Processing，1985，33(3):527-536.