基于協方差矩陣Toeplitz重構的MUSIC算法的波達方位估計方法

王緒虎，侯玉君，金序，田雨，王辛杰，李恩玉

(青島理工大學信息與控制工程學院，山東青島 266520)

0 引言

波達方向(Direction of Arrival,DOA)估計是陣列信號處理的一個重要研究方向[1]。由于水聲信道的復雜性，信號在實際環境傳播中存在多徑傳播和同頻干擾，入射到陣列的信號多為相干信號[2]，傳統的DOA 估計方法性能將會下降，甚至完全失效[3-4]。因此，相干信號的方位估計成為重要的研究方向之一[5]。

解相干方法通常分為兩類，一類是降維處理[6]，如空間平滑方法，這類方法能夠很好地解相干，但是損失了部分陣列的有效孔徑[7]。另一類是非降維處理，如最大似然估計方法[8]，這類方法對相干信號有效，但通常利用多維搜索進行計算，計算量巨大，且求解過程中存在局部最優解的問題，實用性不強。文獻[9]提出了一種基于前向和后向部分托普利茲(Toeplitz)矩陣重構的改進方位估計方法，利用樣本協方差矩陣的半行來重構數據矩陣，解決了因入射信號與角度間存在相位差造成DOA估計性能降低的問題。文獻[10]提出了一種基于多重Toeplitz矩陣集重構的相干信號估計方法，該方法利用陣元接收數據構造托Toeplitz矩陣，通過厄爾米特(Hermitian)轉置矩陣修正和正反向平滑運算得到滿秩矩陣實現解相干，具有較好的估計精度和成功分辨率。文獻[11]提出了一種雙目標相干源窄帶信號波達方向估計方法，該方法不受信號相關性的影響，所需數據量少，提高估計精度的同時降低了計算量。文獻[12]針對基于子空間方法不能解決信號數量大于傳感器數量及基于稀疏性方法總是受到基失配影響的問題，提出了一種協方差矩陣重構方法，提高了DOA估計的速度和精度。文獻[13]提出了一種基于均勻線陣改進的多重信號分類(Modify Multiple Signal Classification,MMUSIC)算法，在不增加計算量的前提下有效降低了信號間的相干系數，提升了MUSIC 算法的估計性能，但是該方法僅適用于線列陣。文獻[14]提出了一種基于重構噪聲子空間(Reconstructed Noise Subspace,RNS)的相干信號估計方法，并且在低信噪比、少快拍下估計精度優于MMUSIC 方法，但是當信號入射間隔較小或者存在多個相干信源時，DOA 估計性能會嚴重降低。

因此，本文首先根據樣本協方差矩陣及其翻轉矩陣構造新協方差矩陣，并利用新協方差矩陣構造Toeplitz 矩陣，然后利用MUSIC 算法進行方位估計。仿真實驗和湖上實驗數據處理結果表明，與MMUSIC 方法和基于RNS 的估計方法相比，本文所提方法提高了DOA估計精度和方位分辨能力。

1 信號模型

假設有K個遠場窄帶信號入射到由M個陣元構成的均勻線列陣上(K

圖1 均勻線陣模型Fig.1 Uniform linear array model

以左側第一個陣元為參考，則此時的陣列導向矢量為

第m個陣元接收到的信號可以表示為

其中：m=1,2,…,M，sk(t)(k=1,2,…,K)是第k個入射到均勻線列陣上的信源；λ是信號波長，nm(t)是高斯白噪聲，θk是第k個信源的入射角度，[ ?]T表示矩陣的轉置。

此時的陣列數據接收矩陣形式可以表示為

其中：

由式(3)中的X(t)，可得陣列輸出信號的協方差矩陣為

其中：RS=E[S(t)SH(t)]表示信號的協方差矩陣，σ2表示噪聲方差，IM表示M階單位矩陣，[ ?]H表示矩陣的共軛轉置。

2 算法原理

2.1 MUSIC算法

對式(8)中的協方差矩陣RX進行特征值分解可得：

其中：US表示由K個大特征值對應的特征向量組成的信號子空間，ΣS表示大特征值對應的對角矩陣，UN表示由M-K個小特征值對應的特征向量組成的噪聲子空間，ΣN表示小特征值對應的對角矩陣。因為陣列流型矩陣A(θ)中的每一個列向量都與噪聲子空間正交，所以理想狀態下有：

對式(11)的空間譜進行譜峰搜索得到DOA的估計值。當入射信號不相干時，該方法具有較好的方位估計性能，當入射信號相干時，該方法性能嚴重降低，不能正確地估計出入射信號的方位。

2.2 MMUSIC算法

由于MUSIC 算法無法對相干信號進行方位估計，研究者在該方法基礎上提出了一種能夠估計相干信號的改進MUSIC 算法[13]。改進算法對陣列數據接收矩陣X(t)預處理，得到一個新的數據接收矩陣Y(t)，表達式為

其中：J是M×M階副對角線元素為1、其余元素均為0的反對角矩陣；X*(t)為陣列接收數據矩陣X(t)共軛矩陣；[ ?]*表示矩陣的共軛。

此時，信號的協方差矩陣可表示為

對協方差矩陣RX、RY取平均得到新的協方差矩陣：

然后對新協方差矩陣R進行特征分解，得到新信號子空間和新噪聲子空間，將新噪聲子空間代入式(11)計算空間譜，然后對空間譜進行譜峰搜索，譜峰對應的角度即為入射信號方位的估計值。該方法可以有效地估計出兩個相干信號的入射方位。

2.3 RNSMUSIC算法

MMUSIC算法只利用了X(t)和Y(t)的自協方差信息，重構噪聲子空間MUSIC(Reconstructed Noise Subspace Multiple Signal Classification,RNSMUSIC)算法[14]不僅利用了X(t)和Y(t)自相關信息，還利用了兩者之間的互協方差信息。

利用X(t)和Y(t)的自協方差矩陣、互協方差矩陣構造增廣矩陣RA，計算方法為

其中：RXY=E[X(t)YH(t)]是X(t)和Y(t)的互協方差矩陣，RYX=E[Y(t)XH(t)]是Y(t)和X(t)的互協方差矩陣。

對RA做奇異值分解，得到特征值矩陣Z和特征向量矩陣U：

對式(14)中的新協方差矩陣R做奇異值分解，得到特征值矩陣Z1和特征向量矩陣U1：

將特征值矩陣Z和Z1的對角線元素按順序比較，選取較大值構成新的特征值矩陣Z0：

其中：λ1≥λ2≥…≥λM。

根據特征值矩陣Z0重構新的特征向量矩陣U0：

從矩陣U0中選擇第個K+1列到第M列構成新的噪聲子空間Un：

然后根據Un構造空間譜函數：

最后，通過譜峰搜索確定DOA的估計值。

雖然MMUSIC 算法和RNSMUSIC 算法都能夠準確地估計兩個相干信號的波達方向，但是相干信號的個數超過兩個時，這兩種方法的方位估計性能下降。

2.4 本文方法

由于噪聲的影響，破壞了協方差矩陣的Hermi‐tian 特性，導致方位估計性能下降，針對該問題，本文通過對接收數據協方差矩陣進行Toeplitz重構來改善方位估計性能。本文方法選取接收數據協方差矩陣的第一行和第一列，構造具有Hermitian 特性的新協方差矩陣，然后利用MUSIC 方法估計目標方位。該方法可以估計出相干信號，并且不降低空間分辨力[15]。

協方差矩陣R中的第m行、n列的元素可以表示為

分別取R中的第1 列和第1 行為參考元素，記為e1，e2：

對式(25)中的空間譜函數進行譜峰搜索，得到DOA估計值。

綜上所述，本文方法的具體實現過程如下：

(1)根據式(3)構造信號協方差矩陣RX。

(2)設置交換矩陣J，根據式(12)構造新的數據接收矩陣Y(t)和協方差矩陣RY。

(3)根據式(14)得到新的協方差矩陣R。

(4)對R進行Toeplitz重構得到。

(6)對式(25)中的空間譜函數進行譜峰搜索得到DOA的估計值。

3 仿真分析

仿真實驗中采用陣元數為8的均勻線陣，陣元間距為0.5 m，兩個中心頻率為1.5 kHz的遠場窄帶信號分別從0°和15°的方向入射到陣列。本文分別采用MMUSIC算法、RNSMUSIC算法、MUSIC算法以及本文方法進行DOA 估計。為了驗證本文方法具有較高的分辨能力，將信噪比固定為0 dB，采樣快拍數分別設置為200和50，圖2、3給出了兩個不同快拍數時的空間譜圖。

由圖2和圖3可以看出，在信噪比為0 dB、快拍數為200時，MUSIC方法無法解相干，不能估計出兩個相干信號，MMUSIC 方法、RNSMUSIC 方法和本文方法都可以成功分辨兩個相干信號，但是本文方法的譜峰更加尖銳，旁瓣更低。在信噪比為0 dB，快拍數為50 時，MMUSIC 方法產生偽峰干擾，估計性能嚴重下降，RNSMUSIC 方法和本文方法可以成功分辨兩個相干信號，并且本文方法估計性能更好。

圖3 快拍數為50的DOA估計空間譜圖Fig.3 Spatial spectrum of DOA estimation with 50 snapshots

將兩個信號的入射角度分別設置為0°和7°，采樣快拍數為200，其他仿真條件不變。圖4 給出了四種方法的空間譜圖。從圖4可以看出，本文方法方位分辨力最好，方位估計精度最高。

圖4 入射信號間隔較小時的DOA估計空間譜圖Fig.4 Spatial spectrum of DOA estimation when the signal interval is small

將信號的入射角度分別設為0°、15°、25°，采樣快拍數為200，其他條件不變，得到圖5所示的空間譜圖。從圖5中可以看出，本文方法可以估計更多數目的相干源。

圖5 信源個數為3時的DOA估計空間譜圖Fig.5 Spatial spectrum of DOA estimation when the number of sources is 3

圖6 三種方法分辨成功概率隨信噪比的變化Fig.6 Variation curves of the success probabilities of three methods with signal to noise ratio

圖7 三種方法分辨成功概率隨快拍數的變化Fig.7 Variation curves of the success probabilities of three methods with the number of snapshots

圖8 三種方法分辨成功概率隨角度間隔的變化Fig.8 Variation curves of the success probabilities of three methods with the angle interval

從圖6可以看出，隨著信噪比的增加，信號的分辨成功概率也增加。在信噪比為3 dB 時，本文方法的分辨成功概率達到100%，與MMUSIC方法和RNSMUSIC 方法相比，可分辨信噪比分別降低了6 dB和8 dB，提高了空間方位分辨能力。

從圖7可以看出，本文方法優于其他方法，在快拍數為50 時，本文方法的分辨成功概率能夠達到86.5%，MMUSIC 方法和RNSMUSIC 方法的分辨成功概率分別為18.5%和40.5%。本文方法在快拍數為300 時，成功概率達到了100%，而MMUSIC 方法和RNSMUSIC 方法，在快拍數達到400時，成功概率依然在90%以下。因此，本文方法具有更高的穩定性，在快拍數較少時依然具有較好的方位估計性能。

從圖8 可以看出，隨著角度間隔的不斷增大，三種方法的分辨成功概率也在不斷增加。在入射角度間隔為8°時，本文方法的分辨成功概率能夠達到100%，MMUSIC方法和RNSMUSIC方法的分辨成功概率分別為95%和97%，在角度間隔為12°時，這兩種方法的分辨概率仍然低于100%。因此，本文方法在角度間隔較小時具有更好的分辨能力。

4 實驗分析

本節對湖上實驗的數據進行處理，驗證本文方法在實際應用環境中的有效性。實驗時，湖面布放兩只實驗船，發射船的頭部和尾部分別懸掛兩個發射換能器，發射同頻率的CW脈沖信號，接收船吊掛一個5 個陣元的均勻水聽器線陣接收聲源信號，陣元間距按發射信號的半波長配置，兩只小船的距離約為300 m，船長6 m，由幾何關系可以分別推算出兩個聲源入射角度約為14°和20°。圖9給出了湖上實驗設備的位置示意圖。

圖9 實驗設備布設圖Fig.9 Layout diagram of test devices

從水聽器采集的湖試數據中截取了20 000個樣本數據，圖10為5個陣元接收的數據結果圖。分別采用本文方法、MMUSIC 方法、RNSMUSIC 方法以及MUSIC 方法對實驗采集到的數據進行處理，得到的方位歷程圖如圖11所示。

圖10 陣元的實際接收數據Fig.10 The actual data received by the array

圖11 四種方法的方位歷程圖Fig.11 Target orientation course charts obtained by four methods

從圖11 中的4 個子圖可以看出，本文方法在DOA 估計時充分利用了陣列接收數據，可以清楚地分辨兩個相干信號。MMUSIC方法在DOA估計時則沒有充分利用陣列接收數據，導致該方法無法成功分辨出兩個相干信號。RNSMUSIC 方法在DOA 估計時雖然充分利用了陣列接收數據，但是由于信號入射角度間隔很小，無法分辨出兩個相干信號。MUSIC 方法無法對相干信號進行解相干，不能有效地估計信號的DOA。從而驗證了本文方法在相干信號DOA估計中的優異性能。

5 結論

本文提出了一種基于協方差矩陣Toeplitz 重構的MUSIC 算法的波達方位估計方法。該方法首先根據陣列接收數據的協方差矩陣及其翻轉矩陣來構造新協方差矩陣，并利用新協方差矩陣構造To‐eplitz 矩陣，使陣列接收數據得以充分利用，然后通過MUSIC 方法進行空間譜估計。仿真和實驗數據處理結果表明，與MMUSIC 方法、RNSMUSIC方法以及MUSIC方法相比，本文方法在低信噪比、少快拍以及入射角度間隔較小條件下的DOA 估計性能更好，并且能夠對3個以上的多個相干信號進行精確的方位估計，在實際水下環境應用場景中也具有較高的方位分辨率和較好的方位估計性能，具有較大的工程應用價值。