周 圍，冉靖萱，彭 洋，陳星宇，馬茂瓊

（1.重慶郵電大學 光電工程學院，重慶 400065；2.重慶郵電大學 移動通信技術重慶市重點實驗室，重慶 400065）

0 引 言

單基地大規(guī)模多輸入多輸出（Multiple?Input Multiple?Output，MIMO）雷達是近年來在通信領域中廣泛使用的新型雷達，其憑借著高速、穩(wěn)定等優(yōu)點，成為了5G 通信系統(tǒng)的核心技術。MIMO 雷達使用多個陣元同步發(fā)射和接收信號以執(zhí)行集中信號處理，這對于在諸如提升信道容量、提高信號檢測和目標參數(shù)估計性能等方面具有顯著的優(yōu)勢。

波達方向（Direction of Arrival，DOA）估計在通信、醫(yī)學成像、天文成像、地震學、聲吶等領域有著廣泛的應用。由于單基地MIMO 雷達具有同時估計相干和非相干目標的能力，因此DOA 估計是MIMO 雷達研究中的一個重要方向。近年來所提出的MIMO 雷達DOA 估計算法大都基于陣元間隔為半波長的均勻線性陣列，在這種情況下總的可分辨目標數(shù)為接收陣元數(shù)減1，可分辨的相干目標為發(fā)射陣元數(shù)，并且在大規(guī)模MIMO 陣列中很難在有限的空間內部署非常多的陣元，因此本文在MIMO 陣列中引入互質陣列，通過研究互質陣列的特點，用更少的陣元達到均勻直線MIMO 陣列需要更多陣元才能估計的目標數(shù)。

互質陣列是一種非均勻稀疏陣列，相比均勻陣列，互質陣列可以對入射信號進行欠采樣，從而突破陣元之間的間距限制，并且陣列可以獲得遠超陣元個數(shù)的自由度（Degree of Freedom，DOF），但是這些方法中并沒有考慮目標相干的情況。文獻[11]對單基地MIMO 雷達進行了研究，得出了可以成功估計的相干目標和非相干目標的數(shù)量與陣元之間的關系。文獻[12]將互質陣列引入MIMO 雷達，結合ESPRIT 算法進行DOA 估計，但該算法只能估計非相干目標。文獻[13]介紹了一種互質MIMO 陣列的幾何結構，文章從互質MIMO 陣列的原始數(shù)據(jù)中模擬出和協(xié)陣列（Sum Co?Array）的協(xié)方差矩陣，從而大大增加了互質陣列的DOF。

為了提高DOA 估計的性能，并且盡可能地降低算法的復雜度，本文重點針對互質MIMO 陣列估計陣列的幾何模型進行了研究，找到了一種和協(xié)陣列，并利用多頻操作對陣元進行擴展，提出基于MUSIC 算法的發(fā)射/接收和協(xié)陣列DOA 估計方法（T/R Co?Array DOA Estimation Method Based on MUSIC，TR?CA?MUSIC），通過重排列數(shù)據(jù)矩陣構造MUSIC 算法譜峰搜索需要虛擬的數(shù)據(jù)矩陣。對于個發(fā)射陣元、2-1 個接收陣元的互質MIMO 陣來說，該和協(xié)陣具有+-1 個連續(xù)陣元和2個中間只有一個孔的半連續(xù)陣元。通過改變發(fā)射陣元的信號發(fā)射頻率，將半連續(xù)陣元中的兩個孔洞用附近的陣元補齊，使得整個連續(xù)陣列陣元數(shù)目擴大為+3-1，該陣列最多可以估計+3-2 個非相干信號、(+3-1) 2 個相干信號。同時，為了避免MUSIC 算法中特征值分解帶來的巨大運算量，本文引入傳播算子，提出基于傳播算子的發(fā)射/接收和協(xié)陣 列DOA 估 計 方 法（T/R Co?Array DOA Estimation Method Based on PM，TR?CA?PM）。其中：TR?CA?MUSIC方法的運算量較大，但可以同時準確估計相干和非相干目標；TR?CA?PM 算法在少目標情況下運算量大大降低，但其相干目標估計能力較弱，比較適合少量非相干目標的DOA 估計。最后本文分別給出了兩種算法的空間譜圖和均方根誤差圖，證明了兩種算法目標數(shù)的估計能力，并驗證了在信噪比（SNR）大于-5 dB 時，TR?CA?PM 的性能與TR?CA?MUSIC 接近。

1 系統(tǒng)模型

1.1 信號模型

假設和互為質數(shù)，且。單基地互質MIMO 陣列如圖1 所示，該陣列包含了個發(fā)射陣元和2-1 個接收陣元，其發(fā)射和接收陣元的位置可表示為：

圖1 單基地互質MIMO 陣列幾何構型

考慮到等效虛擬陣列為均勻線陣，所以定義=2，為發(fā)射陣列發(fā)射的窄帶信號的波長。

假設有個非相關目標和個相關目標位于陣列的遠場，陣列的個發(fā)射器中的每個發(fā)射陣元都發(fā)射一個窄帶信號照射目標，然后通過所有的2-1 個接收器收集所有目標的反射信號，從而形成一個(2-1) ×維的信號矩陣。

式中：第(,)項表示第個發(fā)射器發(fā)射并由第個接收器接收到的信號矢量；=c，為信號的波數(shù)，為發(fā)射陣元的發(fā)射角頻率（信號角頻率），c 為光速；r,t分別為接收陣元和發(fā)射陣元的位置；θ為第個目標的真實角度；ˉ是均值為零、方差為的高斯白噪聲；ˉ可以表示為：

式中：α是信號的振幅；β是非相干目標的相位；β是相干目標的相位。因此的協(xié)方差矩陣為：

式中：[]表示期望運算；H 表示共軛轉置。

1.2 和協(xié)陣列模型

傳統(tǒng)MUSIC 算法只需要利用1.1 節(jié)中的協(xié)方差即可進行DOA 估計，但由于互質MIMO 陣列是稀疏陣列，陣元間距大于半波長，直接進行DOA 估計會產生相位模糊。

根據(jù)式（2）可知，原始數(shù)據(jù)的相位由發(fā)射器位置的和確定，所以可以找到想要的等價虛擬陣列，然后通過重新排列原始數(shù)據(jù)矩陣構造虛擬數(shù)據(jù)矩陣即可估計DOA。

圖1 的和協(xié)陣列如圖2 所示，可以看出該陣列并不是一個完整的連續(xù)線性陣，其中有一些孔洞并沒有填充陣元，但可以將其劃分為三個部分：陣列正中的連續(xù)子陣；緊挨連續(xù)子陣兩旁的半連續(xù)子陣；最外邊的非連續(xù)子陣。連續(xù)子陣占了該和協(xié)陣列的很大一部分，經計算，連續(xù)子陣包含+-1 個陣元，半連續(xù)子陣包含2-2 個陣元，從而，該互質MIMO 陣列利用+數(shù)量級的陣元數(shù)目達到了×數(shù)量級的自由度。

圖2 單基地互質MIMO 陣列的和協(xié)陣列幾何構型

1.3 多頻操作延長連續(xù)子陣

如果將半連續(xù)子陣中的兩個孔填上，則和協(xié)陣列就會多出2個連續(xù)陣元，自由度也將增加2。在單頻率模式下虛擬形成的陣列孔可以通過多頻率操作填充孔洞，以形成更長的連續(xù)子陣。由式（2）可得：

令工作頻率ω=，則當工作頻率為ω時，陣列流形的相位可以表示為：

式中：(r+t)是頻率為ω時和協(xié)陣列的位置，相比初始頻率時的位置r+t，可以將新陣列視為原陣列按比例延伸倍。

圖3 三種工作頻率下的和協(xié)陣列幾何構型

2 DOA 估計

2.1 TR?CA?MUSIC 方法DOA 估計

擴展后連續(xù)子陣的位置可以表示為：

假設從連續(xù)子陣中取個發(fā)射陣元和個接收陣元，則發(fā)射和接收陣元的位置可以表示為：

式中α為任意遞增的正整數(shù)，也就是說相鄰兩個發(fā)射陣元的間距是不固定的。當滿足式（10）時，虛擬陣列和連續(xù)子陣近似等價。

式 中 ： ?∈[1,2,…,]；∈[1,2,…,+3-1]；∈[1,2,…,]。特別地，當α=時，連續(xù)子陣列包含的虛擬陣元數(shù)目最多，所以自由度最大，此時發(fā)射陣元和接收陣元的位置可表示為：

式中為任意正整數(shù)。根據(jù)不同的相干和非相干信號的個數(shù)，可以靈活地設置和的值。

接下來分四步進行DOA 估計：

1）寫出和協(xié)陣列中+3-1 個連續(xù)陣元的數(shù)據(jù)矩陣。

2）構造×維的數(shù)據(jù)矩陣X：

式中X由X的個連續(xù)元素組成。

3）將式（13）代入式（4）得到×維的數(shù)據(jù)協(xié)方差矩陣：

4）對R進行特征值分解，并利用MUSIC 算法進行譜峰搜索。

式 中：E是 協(xié) 方 差 矩 陣R的 噪 聲 子 空 間；() =[1,e ,e ,…,e ]，是個接收陣元的陣列流形向量。

該方法有如下幾個特點：

1）當＞＞1 時，最多可以分辨-1 個目標（包含相干和非相干），其中最多可分辨?zhèn)€相干目標。

2）當＞＞1 時，最多可分辨-1 個任意相干或非相干目標。

3）當=1 時，最多只能分辨-1 個非相干目標，且此時可以分辨的目標數(shù)最多。

2.2 TR?CA?PM 方法DOA 估計

以上的TR?CA?MUSIC 方法涉及到對×維的協(xié)方差矩陣R進行特征值分解，這會導致算法的計算量極大，影響算法的速度。傳播算子算法（Propagation Method，PM）利用矩陣分塊代替特征值分解獲得噪聲分量，可以有效避免特征值分解帶來的巨大運算量，因此本文引入傳播算子算法提出TR?CA?PM 算法。利用傳播算子獲得噪聲子空間的步驟如下：

1）對R進行分塊，有：

式中：[ ? ]表示矩陣左右分塊；為×(+)維矩陣，是R的前+列；為R的剩余項。

2）存在一個(+)×(-(+))維的矩陣，使得=成立，其中就是傳播算子，由最小二乘法得：

4）根據(jù)步驟3）構造譜峰搜索函數(shù)：

2.3 復雜度分析

TR?CA?MUSIC 方法的計算量主要是由協(xié)方差矩陣運算、特征值分解以及譜峰搜索帶來的。假設將乘法運算次數(shù)看作復雜度，則其協(xié)方差運算式（14）的復雜度為(?)，特征值分解復雜度為()，譜峰搜索復雜度為((-(+) )+)( 1800 step )，其中，step 為搜索步長，本文取0.01。兩種方法中的協(xié)方差矩陣運算和譜峰搜索的復雜度相同。所以，TR?CA?PM 方法運算量不同的地方是矩陣運算，其復雜度為((+)+ (+))。經Matlab 計算驗證，當信號數(shù)量+＜34 時，TR?CA?PM 方法復雜度低于TR?CA?MUSIC 方法。特別地，當目標數(shù)量相對陣元數(shù)目越少時，復雜度降低則越明顯，在大規(guī)模MIMO 陣列中，復雜度的降低極為顯著。當目標數(shù)量最多（+=-1）時，兩種方法復雜度相當。

3 仿真結果

3.1 TR?CA?MUSIC 方法空間譜仿真

通過Matlab 2016a 仿真分析TR?CA?MUSIC 方法的相干和非相干目標識別能力。設置互質MIMO 陣列中和分別為3 和4，假設使用全部連續(xù)陣元，目標的信噪比SNR=10 dB，信號的快拍數(shù)Snapshot=1 000。

如圖4 所示，當陣列的發(fā)射陣元數(shù)=10，接收陣元=11 時，方法成功識別了-1=10 個相干目標。圖5顯示了=2，=19時算法的識別效果，其中一共識別了18個目標，包含2個相干目標（左2）和16個非相干目標。如圖6 所示，當=1，=20 時，陣列不能識別相干目標，但此時是識別目標最多的情況，即可以識別-1=19 個非相干目標。

圖4 A=10，B=11 時相干目標的MUSIC 譜圖

圖5 A=2，B=19 時混合目標的MUSIC 譜圖

圖6 A=1，B=20 時非相干目標的MUSIC 譜圖

3.2 TR?CA?PM 方法空間譜仿真

仿真分析TR?CA?PM 方法的相干和非相干信號識別能力，保持，，SNR，Snapshot大小不變。

由圖7 可看出，=1，=20 時TR?CA?PM 算法同樣能成功分辨出最多19 個非相干目標，所以該方法在整體復雜度較低的情況下仍能達到TR?CA?MUSIC 方法的識別效果，但該方法的相干目標分辨能力弱于TR?CA?MUSIC 方法。由圖8 可知，TR?CA?PM 算法僅能識別相干目標的大致角度，因此該方法可用做相干目標角度的快速預處理或者相干目標數(shù)量的快速識別。

圖7 A=1，B=20 時非相干目標的PM 譜圖

圖8 A=2，B=19 時混合目標的PM 譜圖

3.3 均方根誤差分析

本節(jié)的仿真中仿真條件與上節(jié)相同，通過對比TR?CA?MUSIC 方法和TR?CA?PM 方法的均方根誤差（Root Mean Square Error，RMSE）隨信噪比的變化情況來衡量兩種算法估計非相干目標的性能。均方根誤差公式如下：

式中：+為相干和非相干目標個數(shù)之和；為蒙特卡洛實驗次數(shù)；θ為第個信號的第次實驗的DOA 估計值；θ為信號DOA 實際值。

設置蒙特卡洛次數(shù)=1 000，信噪比（SNR）從-11~10 間隔取3，得到兩種算法的仿真圖如圖9 所示?？煽闯霎斝旁氡群涂炫臄?shù)增大時，兩種算法的均方根誤差均隨之減小。TR?CA?PM 算法在SNR≥-5 時均方根誤差非常接近TR?CA?MUSIC 方法，因此TR?CA?PM 算法有著更低的復雜度且具有與TR?CA?MUSIC 相似的準確率。

圖9 A=10，B=11 時兩種方法的RMSE 隨SNR 變化圖

4 結 語

本文通過對互質MIMO 陣列的和協(xié)陣列進行研究，提出了基于MUSIC 的發(fā)射/接收和協(xié)陣列DOA 估計方法。通過附加另外兩種工作頻率，延長了和協(xié)陣列中連續(xù)陣元的數(shù)量，提高了陣列的自由度，用2+-1 個互質MIMO 陣元，實現(xiàn)了最多+3-2 個目標的DOA 估計能力，運用到大規(guī)模MIMO 陣列中時，自由度的提升更加明顯；然后利用數(shù)據(jù)矩陣重構的方式解相干，實現(xiàn)對相干和非相干目標的DOA 估計；最后針對MUSIC 算法中的特征值分解計算量極大的問題，提出了基于PM 的發(fā)射/接收和協(xié)陣列DOA 估計方法。在目標為非相干的情況下，兩種方法均能有效識別目標，仿真證明其性能相當，在識別混合目標時，后者能以估計相干目標的數(shù)量大致估計其角度。此外，兩種方法均可根據(jù)目標的數(shù)量調整發(fā)射/接收和協(xié)陣元數(shù)目，以達到最低的復雜度。