基于插值矩陣技術與二維空間平滑法的近場相干聲源識別研究

單穎春，劉耀光，劉獻棟，何 田

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

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

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

1 聲源信號模型

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

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

當K個聲信號從近場入射時，則:

第(m，n)個傳聲器接收到的信號xmn為:

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

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

式中，

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

當K個聲信號從遠場入射時，傳聲器接收到的信號表示成向量形式為:

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

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

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

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

2.1 插值矩陣技術

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

為了得到近場信號與遠場信號之間的轉換矩陣，分別以間隔Δθ，Δφ對方位角和仰角進行劃分，如圖3所示。

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

式(5)中:

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

對于同樣的信號Svirtual(t)，如果其來自于遠場，則傳聲器陣列接收到的信號為:

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

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

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

對式(2)左乘TH，把實際在近場中接收到的信號Xn(t)轉換為遠場信號Xf(t)，即:

以上步驟為進一步應用空間平滑法去除信號的相關提供了必要的條件。

2.2 二維空間平滑法

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

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

第(m，n)個子陣列接收到的信號為:

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

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

其中，RSS=E［S(t)S(t)H］表示信號的協(xié)方差矩陣，當信號互不相關時是一個對角矩陣，當信號相干時是非對角矩陣。

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

其中:

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

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

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

2.3 聲源信號的功率預測

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

［8］，對應方向向量信號的功率為:

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

3 仿真分析

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

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

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

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

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

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

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

4 結論

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