非均勻線陣的快速求根MUSIC方法?

陳柱學,吳建新

(1.中國電子科技集團公司第三十八研究所，安徽合肥230088;2.西安電子科技大學雷達信號處理國家重點實驗室，陜西西安710071)

0 引言

DOA估計是陣列信號處理的一個重要應用[1-3]。在各種DOA估計方法中，MUSIC方法[4-5]由于它的性能優異而得到了廣泛應用。但傳統MUSIC方法需要進行角度搜索，精細的角度搜索往往導致很大的計算量。為了減少搜索運算，Bhaskar提出了求根 MUSIC方法[6]，它將DOA估計轉化為復多項式求根問題，直接由復根的相位來確定角度。但求根MUSIC方法要求陣列為均勻陣列，對于非均勻陣列，該方法就會失效。在陣元數相同的情況下，非均勻線陣能夠達到比均勻線陣更大的孔徑，提高了角分辨率[7]。因此，非均勻陣列的DOA估計是實際應用中面臨的一個重要問題。Friedlander和Doran針對非均勻陣列分別提出了陣列內插方法[8]和流型分離方法[9]，將非均勻陣列變換成虛擬均勻陣列。于是各種適合于均勻陣列的各種DOA算法，包括ESPRIT方法[10]和求根MUSIC方法都能應用于非均勻陣列。但是陣列內插方法和流型分離方法要求虛擬均勻陣列孔徑大小不小于非均勻陣列孔徑大小的兩倍，否則精度會變差。因此，虛擬均勻陣列一般都比較大，直接采用求根MUSIC方法需要對高階多項式進行求根運算，計算量很大并且數值不穩定。

針對非均勻線陣DOA估計帶來的高計算量問題，本文提出了一種快速的多項式求根MUSIC方法。該方法利用小角度范圍內的導向矢量可以利用低階多項式很好近似的特性，將非均勻線陣DOA估計問題轉化為多組低階實多項式求根問題。由于低階實多項式具有計算量小以及數值穩定性好等優點，而且多組多項式求根可以并行實現，該方法能夠快速地實現非均勻線陣DOA估計，非常適合于工程應用。

1 信號模型

考慮N個陣元的線陣，其第n個陣元相對于第1個陣元的間距為d n，其中d1=0，其陣列示意圖如圖1所示。假設有P個遠場基帶信號入射到該陣列，其中第p個信號s p(k)的入射方向與陣列法線方向夾角為θp。k時刻第n個陣元的接收回波信號x n(k)可以表示為

式中，λ表示波長;w n(k)表示k時刻第n個陣元的高斯白噪聲。將式(1)寫成矢量形式為

圖1 非均勻陣列示意圖

回波數據對應的協方差矩陣可以表示為

一般情況下理想的協方差矩陣是得不到的，需要從數據中估計得到，一般采用如下方式估計:對估計的協方差矩陣進行特征分解，

傳統的M USIC算法[1]構造如下的代價函數:

對于非均勻陣列來說，我們可以通過各種途徑將其變換成均勻陣列。兩種常用的變換方法是陣列插值方法[4]和陣列流型分離方法[5]。這兩種方法可以統一寫為

將式(7)代入到MUSIC倒數譜可以得到

由于aT(θ)是一個等比數列，我們可以直接利用求根MUSIC方法將上述代價函數轉化為復多項式，也就是令并且

上式是關于z的2Q-2階復多項式，對這么高階的復多項式進行求復根運算所需要的計算量很大，并且數值穩定性也不好，很不利于工程實現。比如對于一個長度為10λ的6個陣元的非均勻線陣，其虛擬均勻線陣的半波長陣元數為40，此時式(9)的多項式階數為78階。為了降低非均勻陣列求根MUSIC方法的計算量，需要采用有效的近似方法。

2 快速實多項式求根MUSIC方法

式中，g∈[0，ε]。對應地，落在第k個區間內的頻率對應的導向矢量可以表示為

式中，u k=[e-j2π(M-1)εk… 1 … ej2π(M-1)εk]T，c k(g)表示第k個區間的導向矢量。將式(14)代入到式(10)可以得到第k個區間內的M USIC倒數譜為

式中，7表示Hadamard積;Φ=diag(b)表示對角元素為b的對角矩陣。

對于落在小區間內的導向矢量來說，由于它具有低秩特性，可以采用低階多項式來近似。低階多項式近似需要求解如下的最優化問題:

式中，‖·‖2表示2范數;為L階多項式矢量;T=[t1t2…t L+1]∈C M×(L+1)為多項式擬合矩陣。在式(16)中，多項式階數L影響著近似的精度，L越大，近似精度越高，但后續的多項式求根的計算量也就越大。圖2給出了不同階數情況下的擬合誤差情況。從圖中可以看出，當L大于等于5的情況下，各個頻率的近似誤差可以在-40 dB以下，也就是近似的導向矢量與實際的導向矢量的相關系數可以達到0.999 9。因此，L等于5時就可以獲得很高的近似精度。

圖2 不同階數情況下的擬合誤差(M=16)

式(16)是一個經典的最小二乘問題，其最優解[13-14]為

式中，C=[c(g1)c(g2) …c(g I)]，G=[g1g2…g I]，I為區間內的頻率網格點數。這樣就可以得到低階近似后的導向矢量為

將式(18)代入到式(15)可以得到用低階多項式表示的MUSIC倒數譜

1.2.4 疼痛評分 采用11點數字評分法，以無痛為“0”，最劇烈疼痛為“10”，共11個點描述疼痛程度，分值越高，疼痛程度越重?；颊吒鶕约旱奶弁锤杏X在相應的分值處劃對號。

式中 ，z=[z11…z1N z21…z NN]T，w kl=[w kl11…w kl1N w kl21…w klNN]T。根據式(21)，第k個區間的多項式系數矢量可以寫成

式中，W k=[w k1w k2…w k(L+1)]∈C N2×(L+1)。

最后，令低階多項式式(19)等于零，也就是

式(23)的共軛復根的虛部大小可以作為是否是信源的判斷標準，虛部值越接近于0，表示該共軛復根越有可能是信源對應的共軛復根，而實部就是估計的角度。

為了使這個算法更加清楚，我們給出了該算法的處理流程:

1)離線計算獨立于數據的降維矩陣{W1，W2，…，W K};

2)根據噪聲子空間計算矢量z;

3)對第k個角度區間，根據式(22)計算L+1個多項式系數γk1，γk2，…，γk(L+1);

4)對第k個角度區間的多項式式(23)進行多項式求根;

5)重復步驟3和4直到所有角度區間處理完畢，得到所有的共軛復根，然后根據復根的虛部大小判斷信源對應的根，也就是挑選虛部值最小的P個根，根據根的實部直接計算信源的角度。

根據處理流程，本方法總的計算量為O(4N2P+2(L+1)N2+4ML2)實數乘法，而標準的求根M USIC方法的計算量為O(4N2P+32M3)。由于M遠大于N和L，標準的求根M USIC方法的計算量要遠大于本文方法的計算量。圖3給出了本方法的流程圖。

3 仿真實驗

圖3 本文方法的處理流程圖

圖4 本方法獲得的MUSIC譜對應的共軛復根分布

圖5 均方根誤差隨信噪比的變化關系

圖6給出了標準求根MUSIC方法與本方法的均方根誤差隨快拍數的變化關系，其中本方法采用的多項式階數為5階，信噪比為5dB。從圖中可以看出，兩種方法的均方根誤差隨著快拍數的增加而減小，而且對于不同的快拍數，本方法都可以獲得與標準求根MUSIC方法基本相同的性能。另外，隨著快拍數增加到一定程度，均方根誤差性能的改善逐漸趨緩，這主要是由于求根MUSIC方法的收斂性能是非線性的，快拍數比較少時收斂性能改善快，快拍數比較多時基本收斂，此時的收斂性能改善慢。

圖6 均方根誤差隨快拍數的變化關系

4 結束語

針對非均勻線陣DOA估計帶來的高計算量問題，本文提出了一種快速求根MUSIC方法。該方法利用小角度范圍內的導向矢量可以利用低階多項式很好近似的特性，將非均勻線陣DOA估計問題轉化為多組低階實多項式求根問題。該方法能夠在性能基本保持不變的情況下大大降低標準求根MUSIC算法的計算復雜度，非常適合于工程應用。

[1]王永良，陳輝，彭應寧.空間譜估計理論與算法[M].北京:清華大學出版社，2004.

[2]劉周，李軍，劉紅明，等.MIMO雷達非嚴格正交信號單脈沖測角方法[J].雷達科學與技術，2013，11(3):262-266.LIU Zhou，LI Jun，LIU Hong-ming，et al.Monopulse Angle Measurement of Non-Strictly Orthogonal Signals for MIMO Radar[J].Radar Science and Technology，2013，11(3):262-266.(in Chinese)

[3]林文鳳，周曉青，甘露，等.一種可自動配對的二維波達方向估計算法[J].雷達科學與技術，2013，11(1):33-39.LIN Wen-feng，ZHOU Xiao-qing，GAN Lu，et al.An Automatic Pairing Method for 2-D Direction of Arrival Estimation[J].Radar Science and Technology，2013，11(1):33-39.(in Chinese)

[4]Schmidt R O.Multiple Emitter Location and Signal Parameter Estimation[J].IEEE Trans on Antennas and Propagation，1986，34(3):276-280.

[5]黃秋欽，余嘉，梁炎夏，等.基于改進PASTd的MUSIC算法的DSP實現[J].雷達科學與技術，2010，8(2):183-187.HUANG Qiu-qin，YU Jia，LIANG Yan-xia，et al.The Implementation of Modified PASTd Algorithm Based on DSP System[J].Radar Science and Technology，2010，8(2):183-187.(in Chinese)

[6]RAO B D，HARI K V S.Performance Analysis of Root-Music[J].IEEE Trans on Acoustics，Speech and Signal Processing，1989，37(12):1939-1949.

[7]KASSIS C E，PICHERAL J，MOKBEL C.Advantages of Nonuniform Arrays Using Root-MUSIC[J].Signal Processing，2010，90(2):689-695.

[8]FRIEDLANDER B.The Root-MUSIC Algorithm for Direction Finding with Interpolated Arrays[J].Signal Processing，1993，30(1):15-29.

[9]DORAN M A，DORON M A，WEISS A J.Coherent Wide-Band Processing for Arbitrary Array Geometry[J].IEEE Trans on Signal Processing，1993，41(1):414-417.

[10]GERSHMAN A B，RUBSAMEN M，PESAVENTO M.One-and Two-Dimensional Direction-of-Arrival Estimation:An Overview of Search-Free Techniques[J].Signal Processing，2010，90(5):1338-1349.

[11]WU J X，WANG T，BAO Z.Angle Estimation for Adaptive Linear Array Using PCA-GS-ML Estimator[J].IEEE Trans on Aerospace and Electronic Systems，2013，49(1):670-677.

[12]WU J X，WANG T，BAO Z.Fast Realization of Maximum Likelihood angle Estimation with Small Adaptive Uniform Linear Array[J].IEEE Trans on Antennas and Propagation，2010，58(12):3951-3960.

[13]康兆敏.數值計算方法[M].北京:清華大學出版社，2013.

[14]蘇峰，王昌海，徐征.基于最小二乘的時差定位算法[J].雷達科學與技術，2013，11(6):621-625.SU Feng，WANG Chang-hai，XU Zheng.TDOA Location Algorithms Based on the Least Squares[J].Radar Science and Technology，2013，11(6):621-625.(in Chinese)