劉秋雨，蔣彥雯，*，范紅旗，連紅飛

1.國防科技大學 自動目標識別全國重點實驗室，長沙 410073

2.國防科技大學 電子科學學院，長沙 410073

陣列雷達［1-2］在空間高分辨、參數(shù)估計等方面具有良好的性能，是近年的研究熱點，而波達方向（Direction of Arrival，DOA）估計是陣列雷達一個重要研究方向［3-4］，特別是突破瑞利限的超分辨譜估計方法，例如MUSIC［5］方法、最大似然法［6］（Maximum-Likelihood，ML）、子空間擬合法［7］（Subspace Fitting，SF）、旋轉(zhuǎn)子空間不變法［8］（Estimation of Signal Parameters via Rotational Invariance Techniques，ESPRIT）等，其中，MUSIC 方法因具有角度高分辨、不依賴陣列配置、便于實現(xiàn)等優(yōu)勢備受關(guān)注［9］，然而其高分辨和高精度估計涉及較大的矩陣運算，制約了MUSIC 方法的大規(guī)模普及應(yīng)用［10］。

MUSIC 方法的計算量主要來自子空間分解和譜峰搜索2 個部分［11］，與陣元數(shù)、快拍數(shù)、譜峰搜索次數(shù)成正比關(guān)系，實際系統(tǒng)中可通過增大上述3 個參數(shù)的數(shù)值來提升算法性能，但也提高了系統(tǒng)算力、功耗、硬件成本等要求［12］，降低了算法實時性。

為解決MUSIC 方法復雜度高、計算量大的難題，學者們研究了一系列的改進方法。在子空間分解部分，文獻［13］利用多級維納濾波快速求解信號子空間，但要求已知信號波形；文獻［14］提到利用波束轉(zhuǎn)換矩陣將陣元域輸出數(shù)據(jù)轉(zhuǎn)換到波束域，可以簡化子空間計算，但會丟失部分信息；還有一些無需特征分解的子空間方法，文獻［15-16］分別對互相關(guān)矩陣和傳播算子進行線性運算快速得到噪聲子空間，但估計精度稍差于MUSIC 方法。在譜峰搜索部分，ESPRIT［17］和Root-MUSIC［18］雖避免了譜峰搜索，可直接求解角度，但估計精度不高，性能受限；文獻［19］提出壓縮MUSIC（Compressed MUSIC，C-MUSIC）方法，利用類噪聲子空間簇的交集縮小搜索范圍，但要求陣元數(shù)大于目標數(shù)的2 倍?？傮w而言，目前改進的MUSIC 方法雖然可以降低復雜度，但都是在估計性能和低復雜度間折中進行權(quán)衡取舍。因此，如何緩解MUSIC 估計性能與系統(tǒng)算力、成本、效率之間的矛盾，在算力受限情況下提升DOA 估計性能，是陣列雷達系統(tǒng)研究與應(yīng)用中亟需解決的難點問題。

本文提出集成先驗的低復雜度Nystr?m-MUSIC 方法。首先，Nystr?m 方法［20］是一種利用矩陣低秩性，通過計算部分樣本組成的子協(xié)方差矩陣的特征值分解來估計整個協(xié)方差矩陣的分解，從而獲取近似信號子空間的方法，可應(yīng)用于大數(shù)據(jù)量下的陣列DOA 快速實現(xiàn)。另外，隨著多源傳感器［21］、協(xié)同探測［22］等技術(shù)的發(fā)展，現(xiàn)有系統(tǒng)中存在且包含未被充分利用的目標先驗信息［23］。因此，本文提出集成先驗的低復雜度Nystr?m-MUSIC 方法，一方面，基于Nystr?m 近似估計信號子空間，減小子空間分解的矩陣維度，另一方面，充分利用已有目標角度先驗信息，對譜峰搜索過程進行空間約束，減少譜峰搜索的次數(shù)，從而有效降低了MUSIC 方法復雜度、提高方法運算效率，為后續(xù)方法在實際陣列雷達系統(tǒng)中的多層次推廣應(yīng)用奠定基礎(chǔ)。

全文框架如下：首先，本文在陣列信號模型基礎(chǔ)上，詳細介紹基于Nystr?m 方法的子空間分解和集成先驗的譜峰搜索方法的基本原理；其次，給出了集成先驗的低復雜度Nystr?m-MUSIC 方法處理流程并對比分析了不同方法的理論復雜度；最后，開展不同條件下的仿真實驗，分析驗證本文所提方法的有效性和估計性能。

1 陣列信號模型

圖1 為M個滿足遠場條件的各向同性陣元所組成的均勻直線陣列模型，陣元間距為d，假設(shè)有K個不相干的遠場窄帶信號s(n)=[s1(n)，s2(n)，…，sK(n)]T分別從不同的角度θ=[θ1，θ2，…，θK]T入射到均勻直線陣列，定義波達方向θk，k=1，2，…，K為第k個點目標與陣列法線的夾角。

圖1 均勻直線陣列結(jié)構(gòu)Fig.1 Uniform linear array structure

第m個陣元在第n個快拍的接收數(shù)據(jù)為

式中：n=1，2，…，N，N為快拍數(shù)；fc為信號載頻；sk(n)為第k個入射信號；wm(n)為零均值高斯白噪聲；選定陣元1 為陣列相位中心，則τm，k為第k個信號源入射到陣元m處相對于陣元1 處的時延，可表示為

其中：c為光速。式（1）可改寫為

式中：λ為入射信號波長。

陣元m在N個快拍下的接收數(shù)據(jù)矢量可表示為

則整個陣列在N個快拍下接收數(shù)據(jù)的矩陣形式可表示為

其中：A=[a(θ1)，a(θ2)，…，a(θK)] 為M×K維陣列流形矩陣；為M維導向矢量；S=[s1，s2，…，sK]T為K×N維信號矩陣；sk=[sk(1)，sk(2)，…，sk(N) ]為N維入射信號矢量；W=[w1，w2，…，wM]T為M×N維噪聲矩陣；wm=[wm(1)，wm(2)，…，wm(N) ]是N維零均值加性高斯白噪聲矢量。

2 方法介紹

傳統(tǒng)MUSIC 方法包括子空間分解和譜峰搜索2 個步驟，都需要比較大的計算量。本文在系統(tǒng)算力有限的應(yīng)用約束下，提出了基于Nystr?m方法的子空間分解和集成先驗的譜峰搜索方法，利用Nystr?m 方法低維近似式（5）所示陣列接收數(shù)據(jù)矩陣X=[X1，X2，…，XM]T的協(xié)方差矩陣，并引入先驗信息縮小譜峰搜索范圍，分別降低了子空間分解和譜峰搜索部分的復雜度。因此，本文方法相較于傳統(tǒng)MUSIC 的優(yōu)勢如下：

1）采用Nystr?m 方法，隨機抽取部分陣元的接收數(shù)據(jù)，構(gòu)造低維抽樣矩陣，用其低維協(xié)方差矩陣逼近X高維協(xié)方差矩陣，并代替高維矩陣譜分解，從而降低子空間分解的計算復雜度。

2）相比傳統(tǒng)MUSIC 在雷達大視場大范圍遍歷下的譜峰搜索過程，本文通過集成先驗信息，提前鎖定目標所在區(qū)域，將有效縮小解空間，提高搜索精準性和效率。

2.1 基于Nystr?m 方法的子空間分解

子空間分解中必須計算和分解接收數(shù)據(jù)的協(xié)方差矩陣，才能獲得信號子空間。然而隨著陣元數(shù)的增多，陣列接收數(shù)據(jù)矩陣維度也會增大，此時傳統(tǒng)MUSIC 方法需要大量的矩陣運算。Nystr?m 方法源于積分方程理論［24］，隨后在機器學習領(lǐng)域被用于加速矩陣的分解［20］、特征向量近似［25］等處理。本文將Nystr?m 方法用于MUSIC方法中近似求解協(xié)方差矩陣的特征向量，估計信號子空間，從而降低子空間分解的復雜度。如圖2 所示，基于Nystr?m 方法的子空間分解利用3 個低維矩陣近似高維接收信號的協(xié)方差矩陣，陣列接收數(shù)據(jù)X的協(xié)方差矩陣RXX的Nystr?m 低維近似矩陣可表示為，其中，RXX為M×M維的矩陣；Y為接收數(shù)據(jù)X的沿陣列維隨機抽取得到的抽樣矩陣，大小為Me×N；RXY為M×Me維互協(xié)方差矩陣；為Y的協(xié)方差矩陣RYY的Moore-Penrose 偽逆，大小為Me×Me。

圖2 基于Nystr?m 的矩陣近似Fig.2 Matrix approximation based on Nystr?m

接下來，詳細介紹基于Nystr?m 方法的子空間分解。對于陣列接收數(shù)據(jù)X=[X1，X2，…，XM]T的協(xié)方差矩陣

根據(jù)Nystr?m 方法，可用一組獨立同分布的采樣矢量Y=[y1，y2，…，yMe]T近似獲得ui

因此，利用各陣元間的獨立同分布特性，從陣元{1，2，…，M}中隨機抽取Me個樣本，K≤Me≤M，這Me個陣元的接收數(shù)據(jù)可排列組成Me×N維的低維抽樣矩陣

然后計算Y的協(xié)方差矩陣以及矩陣X與Y的互協(xié)方差矩陣，有

再對RYY進行特征分解

式中：i=1，2，…，Me。

根據(jù)式（7）計算RXX的近似特征向量，即

其歸一化特征向量為

接下來，信號子空間的近似形式可表示為

根據(jù)式（6）～式（14）的推導求解可知，本文的子空間分解方法共需要++(N+K)MMe次復數(shù)乘法，而傳統(tǒng)MUSIC 方法中的子空間分解過程共需要NM2+M3次復數(shù)乘法［26］。由于Me＜M，因此Nystr?m 方法將協(xié)方差計算部分的乘法次數(shù)從NM2降至NMMe，協(xié)方差特征分解部分的乘法次數(shù)從M3降至，雖然本文方法用部分特征向量近似信號子空間額外還需(N+K)MMe次乘法，但兩者綜合比較發(fā)現(xiàn)，本文方法相較于MUSIC 方法復雜度降低。

2.2 集成先驗的譜峰搜索方法

傳統(tǒng)MUSIC 根據(jù)信號子空間跟噪聲子空間的正交性，以一定的步長通過譜峰搜索得到空間各個方向的譜值，將最大空間譜對應(yīng)的方向確定為目標的方向，其中對應(yīng)的偽空間譜可表示為

式中：Uw表示噪聲子空間。最大空間譜對應(yīng)方向即目標方向

圖3 為本文所提出的方法與MUSIC 方法對空間目標進行譜峰搜索的對比圖，圖中藍色虛線部分為傳統(tǒng)MUSIC 搜索空間，搜索范圍為雷達的視場（Field of View，F(xiàn)OV），記作θFOV，解空間為FOV 離散化處理后的方向集合，假設(shè)搜索步長為ΔθFOV，那么此時搜索次數(shù)為

圖3 本文方法與MUSIC 方法譜峰搜索對比圖Fig.3 Comparison chart of spectral peak search with proposed method and MUSIC method

上述MUSIC 方法的譜峰搜索是在雷達視場FOV 內(nèi)實現(xiàn)，該視場定義了雷達在水平方向的角度覆蓋范圍，一般FOV 越大越好，這樣雷達探測范圍就越大。然而大視場下的譜峰搜索范圍廣，搜索次數(shù)增多，將導致譜峰搜索計算量大；而通過縮小搜索空間或增大搜索步長的方法來減少計算量，都可能導致估計性能的下降，其中，縮小搜索空間會造成目標丟失等現(xiàn)象，增大搜索步長會降低MUSIC 的估計精度。

本文方法搜索空間如圖3 中紅色區(qū)域所示，集成已有的目標先驗信息，自適應(yīng)地初始化目標先驗角度θ0，然后將其鄰域U(θ0，qΔθ)作為目標的先驗區(qū)域θprior，其中Δθ為目標先驗角度的估計精度，q為加權(quán)系數(shù)，表示搜索區(qū)域的大小。最大空間譜對應(yīng)方向即目標方向可改寫為

本文方法的搜索范圍為先驗區(qū)域θprior，解空間為先驗區(qū)域離散化處理后的方向集合，假設(shè)搜索步長與MUSIC 一致均為ΔθFOV，則搜索次數(shù)為

鑒于在實際系統(tǒng)中，先驗信息可來源于不同傳感器不同的估計方法，此時算法的復雜度需依具體系統(tǒng)而定，但本著低復雜度的原則，一般選擇計算量少且對整體復雜度影響不大的算法，因此，本文主要對獲取先驗之后的算法復雜度和效率進行對比分析。此時直接對比式（16）和式（17）可知，MUSIC 方法譜峰搜索部分的計算量為JFOV(M+1)(M-K)［27］，本文方法計算量為KM2+JpriorM(M+1)。集成先驗下的θprior?θFOV，則Jprior＜JFOV，故本文方法計算量明顯低于MUSIC 方法的計算量，且本文方法搜索目標的效率在一定程度上不受雷達視場的影響，對大視場中的目標也能快速完成高精度搜索，解決搜索復雜度與性能之間的矛盾。

2.3 方法實現(xiàn)與復雜度分析

集成先驗的低復雜Nystr?m-MUSIC 方法具體實現(xiàn)步驟如算法1 所示。

綜上分析，本文通過降維和縮小解空間的方式降低目標角度估計過程的算法運算量。本文所提出的方法在子空間分解部分需要NMe2+Me3+(N+K)MMe次復數(shù)乘法，譜峰搜索部分需要KM2+JpriorM(M+1)次復數(shù)乘法，所以本文方法總復雜度為

與MUSIC 方法復雜度對比如表1 所示。表1第2～4 行為理論復雜度對比行，后4 行在給定參數(shù)下定量對比了2 種方法復雜度，此假設(shè)條件下，MUSIC 方法約需要4.8×106次復數(shù)乘法，本文方法約需要5.6×105次復數(shù)乘法，整體相較于MUSIC 方法復雜度降低了8 倍，其中子空間分解部分降低約5 倍，譜峰搜索部分在集成先驗下可降低約一個數(shù)量級。

表1 計算復雜度對比Table 1 Comparison of computational complexity

3 仿真分析

本節(jié)將對本文所提出的集成先驗的低復雜度Nystr?m-MUSIC 方法和傳統(tǒng)MUSIC 方法進行仿真對比分析，分別設(shè)置多目標DOA 估計仿真、不同信噪比條件下DOA 估計性能仿真以及方法運行時長仿真3 個仿真實驗。

3.1 多目標DOA 估計仿真

本節(jié)主要評估驗證本文方法對多目標DOA的估計能力。實驗參數(shù)設(shè)置如表2 所示。

表2 多目標DOA 估計參數(shù)設(shè)置Table 2 Parameter settings of multi-target DOA estimation

圖4 為本文所提方法在不同抽樣陣元數(shù)情況下和MUSIC 方法、512 點FFT 方法對目標角度的估計結(jié)果。圖4 中FFT 方法歸一化幅度峰值對應(yīng)的角度為0°、10.69°，可見FFT 方法無法分辨2 個靠得很近的目標；而讀取本文方法和MUSIC方法的歸一化幅度峰值，對應(yīng)的角度分別為0°、10°、11.5°，即為2 種方法對目標角度的估計值，與表2 中理論值一致，表明本文方法與MUSIC方法均能實現(xiàn)對目標的高分辨角度估計，驗證了所提方法具備優(yōu)良的估計目標的能力。

圖4 目標角度估計結(jié)果Fig.4 Results of target angle estimation

3.2 不同信噪比條件下DOA 估計性能分析

根據(jù)2.2 節(jié)的分析，本文方法和MUSIC 方法的DOA 估計性能與譜峰搜索空間范圍及搜索步長密切相關(guān)。本節(jié)在圖3 搜索示意圖的基礎(chǔ)上，根據(jù)先驗信息精度的高低，進一步設(shè)置了2 個實驗分別對比分析DOA 估計性能隨信噪比的變化關(guān)系，信噪比變化區(qū)間設(shè)置在-10～15 dB，其中，實驗1 為同等搜索次數(shù)的條件下，評估本文所提方法與MUSIC 方法對目標DOA 的估計性能，實驗2 則在同等搜索步長的條件下分析2 種方法的估計性能。

本文采用均方根誤差（Root Mean Square Error，RMSE）指標來對估計性能進行定量評估。RMSE 定義為

式中：L表示蒙特卡洛實驗次數(shù)；為第l次實驗中第k個信號θk的估計值。

3.2.1 實驗1：同等搜索次數(shù)下的估計性能

在2.2 節(jié)中，搜索步長離散劃分空間角度會影響譜峰搜索的精度，為分析2 種方法對目標的估計性能，本節(jié)將固定譜峰搜索次數(shù)JFOV=Jprior=101，當先驗為同一傳感器的快速算法，如FFT 等提供的低精度角度信息時，可在同等搜索次數(shù)下采用更小的搜索步長對先驗區(qū)域進行求解。譜峰搜索示意圖如圖5 所示，表3 為相應(yīng)的信噪比及搜索空間參數(shù)，其他陣列系統(tǒng)參數(shù)如表2 所示。

表3 實驗1 參數(shù)設(shè)置Table 3 Parameter settings of Experiment 1

圖5 同等搜索次數(shù)下的譜峰搜索示意圖Fig.5 Diagram of spectral peak search with the same search times

圖6 為同等搜索次數(shù)條件下，經(jīng)過100 次蒙特卡洛實驗，MUSIC 方法和4 種不同抽取陣元數(shù)下本文所提方法對目標方位的估計精度結(jié)果。從圖6 中可以看出，隨著信噪比的增大，估計誤差逐漸減小，最后很快在高信噪比時收斂，其中，本文方法估計精度在0 dB 時收斂到0.03°，MUSIC方法估計誤差在0.29°左右，由此發(fā)現(xiàn)，本文方法估計精度比同等實驗條件下的MUSIC 方法提高了近10 倍。因此，通過引入先驗角度信息縮小空間搜索范圍，在同等搜索次數(shù)下，本文方法譜峰搜索更加緊密，避免了因搜索步長離散劃分空間角度帶來的失配問題對估計結(jié)果產(chǎn)生較大誤差，實現(xiàn)了比傳統(tǒng)MUSIC 方法更高的估計精度。

3.2.2 實驗2：同等搜索步長下的估計性能

由2.3 節(jié)方法復雜度分析可知，算法復雜度與譜峰搜索的次數(shù)成正比，因此，本節(jié)主要目的是在相同搜索精度的條件下，分析2 種方法在不同信噪比下對目標的估計性能。當先驗是由光學傳感器等提供的高精度角度信息時，則需在此精度上進一步縮小譜峰搜索步長，本次實驗設(shè)置傳統(tǒng)MUSIC 方法和本文所提方法的搜索步長均為0.01°，搜索次數(shù)分別為JFOV=10 001、Jprior=1 001，信噪比、搜索范圍等其他參數(shù)設(shè)置如表3，同等搜索步長下譜峰搜索示意如圖7 所示。

圖7 同等搜索步長下的譜峰搜索示意圖Fig.7 Diagram of spectral peak search with the same search step

圖8 為同等搜索步長下，2 種方法對目標方位的估計精度結(jié)果。總體來看，如圖8 所示的曲線，2 種方法隨著信噪比的增大，估計誤差逐漸減小，表明角度估計精度隨信噪比增加在不斷提高。雖然，本文所提方法隨著抽樣陣元數(shù)的減少，估計誤差在低信噪比下逐步增加，但最大增幅僅為0.013°，性能下降??；而且，當信噪比上升至0 dB之后，所提方法的RMSE 與傳統(tǒng)MUSIC 方法的RMSE 曲線幾乎重合，均接近理論極限CRLB，低至10-3量級。因此，在同等搜索步長下，本方法利用16 個抽樣陣元和局部先驗區(qū)域，就實現(xiàn)與64 個陣元大視場下MUSIC 方法相當?shù)母呔裙烙嫛?/p>

圖8 實驗2 估計精度結(jié)果Fig.8 Estimation accuracy result of Experiment 2

結(jié)合2.3 節(jié)中表1 所示理論復雜度，在此實驗條件下定量分析2 種方法的復雜度，MUSIC 方法約需要4.2×107次復數(shù)乘法，本文方法約需要4.3×106次復數(shù)乘法。因此，綜合上述實驗分析，驗證了本文提出的方法能以更低的復雜度實現(xiàn)對目標的高精度估計。

3.3 方法運行時長分析

進一步地，本文引入方法運行時間降低度Ts，表示所提方法相對于MUSIC 方法平均運行時長減少的百分比，其定義為

式中：Tproposed為本文所提出的方法在一定信噪比下進行100 次蒙特卡洛實驗的平均運行時長；TMUSIC為MUSIC 方法在相同條件下的平均運行時長。

本節(jié)實驗主要分析對比本方法和MUSIC 方法在估計位于任意位置的單個點目標角度所用時長，同時為供實際系統(tǒng)算法參考提供一定先驗方法和實驗依據(jù)，本節(jié)實驗引入FFT 說明獲取先驗對算法整體效率的影響，實驗設(shè)置隨機抽取陣元數(shù)依次為1～64 變化，MUSIC 方法和本文所提方法的搜索步長均為0.1°，搜索范圍分別為-50°～50° 和0°～10°，搜索次數(shù)分別為JFOV=1 001、Jprior=101，F(xiàn)FT 點數(shù)為256，其他實驗參數(shù)同表2。經(jīng)過100 次蒙特卡洛實驗，得到2 種方法平均運行時長如表4 所示，不同抽取陣元數(shù)時的運行時間降低度如圖9 所示。

表4 2 種方法平均運行時長和運行時間降低度Table 4 Average running time and running time reduction degree of the two methods

圖9 運行時間降低度隨抽取陣元數(shù)的變化曲線Fig.9 Curves of variation of running time reduction degree with the number of extracted array elements

圖9 中，隨著抽取陣元數(shù)的增多，本文方法運行時間越長。當不考慮獲取先驗的時間，本文方法運行時間降低度從84.9% 逐漸降至80.7%；當考慮先驗獲取的時間，方法運行時間降低度從81.4% 逐漸降至79%。結(jié)合表5 整體而言，二者隨抽取陣元數(shù)下降趨勢一致，獲取先驗的過程使運行時間增長約0.3 ms，運行時間降低度減少2%左右，但始終保持著減少80%左右的運行時間，統(tǒng)計結(jié)果表明，在考慮和不考慮獲取先驗情況下的運行時長和運行時間降低度均相差不大并處于同一量級，但其提供的角度先驗卻可大大降低搜索時間。因此，本文方法通過降低算法復雜度減少了算法運行時長，顯著提高了目標DOA 估計的效率。

4 結(jié)論

1）針對陣列雷達系統(tǒng)中MUSIC 方法估計性能和算法復雜度之間相互制約的問題，本文提出一種集成先驗的低復雜度Nystr?m-MUSIC方法。所提方法基于Nystr?m 近似對子空間分解矩陣進行降維處理，并引入先驗角度信息縮小譜峰搜索范圍，有效降低了MUSIC 方法的復雜度。

2）通過開展多目標DOA 估計和不同信噪比下DOA 估計對比，驗證了本文所提方法的有效性，并發(fā)現(xiàn)在具備與傳統(tǒng)MUSIC 方法同等估計精度的條件下，當抽取陣元數(shù)僅為總數(shù)的25%時，本方法的復雜度降低約8 倍、運行時長降低80%以上，能顯著提高目標DOA 估計效率。

3）綜合本文理論分析和仿真實驗，在集成不同精度的先驗信息下，本文所提方法在具備高估計性能和低復雜度的同時，能快速實現(xiàn)對目標方位的高精度估計，對低算力約束下的陣列雷達系統(tǒng)性能提升具有重要意義。