陳柱學(xué),吳建新

(1.中國(guó)電子科技集團(tuán)公司第三十八研究所，安徽合肥230088;2.西安電子科技大學(xué)雷達(dá)信號(hào)處理國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西西安710071)

0 引言

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

針對(duì)非均勻線陣DOA估計(jì)帶來的高計(jì)算量問題，本文提出了一種快速的多項(xiàng)式求根MUSIC方法。該方法利用小角度范圍內(nèi)的導(dǎo)向矢量可以利用低階多項(xiàng)式很好近似的特性，將非均勻線陣DOA估計(jì)問題轉(zhuǎn)化為多組低階實(shí)多項(xiàng)式求根問題。由于低階實(shí)多項(xiàng)式具有計(jì)算量小以及數(shù)值穩(wěn)定性好等優(yōu)點(diǎn)，而且多組多項(xiàng)式求根可以并行實(shí)現(xiàn)，該方法能夠快速地實(shí)現(xiàn)非均勻線陣DOA估計(jì)，非常適合于工程應(yīng)用。

1 信號(hào)模型

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

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

圖1 非均勻陣列示意圖

回波數(shù)據(jù)對(duì)應(yīng)的協(xié)方差矩陣可以表示為

一般情況下理想的協(xié)方差矩陣是得不到的，需要從數(shù)據(jù)中估計(jì)得到，一般采用如下方式估計(jì):對(duì)估計(jì)的協(xié)方差矩陣進(jìn)行特征分解，

傳統(tǒng)的M USIC算法[1]構(gòu)造如下的代價(jià)函數(shù):

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

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

由于aT(θ)是一個(gè)等比數(shù)列，我們可以直接利用求根MUSIC方法將上述代價(jià)函數(shù)轉(zhuǎn)化為復(fù)多項(xiàng)式，也就是令并且

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

2 快速實(shí)多項(xiàng)式求根MUSIC方法

式中，g∈[0，ε]。對(duì)應(yīng)地，落在第k個(gè)區(qū)間內(nèi)的頻率對(duì)應(yīng)的導(dǎo)向矢量可以表示為

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

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

對(duì)于落在小區(qū)間內(nèi)的導(dǎo)向矢量來說，由于它具有低秩特性，可以采用低階多項(xiàng)式來近似。低階多項(xiàng)式近似需要求解如下的最優(yōu)化問題:

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

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

式(16)是一個(gè)經(jīng)典的最小二乘問題，其最優(yōu)解[13-14]為

式中，C=[c(g1)c(g2) …c(g I)]，G=[g1g2…g I]，I為區(qū)間內(nèi)的頻率網(wǎng)格點(diǎn)數(shù)。這樣就可以得到低階近似后的導(dǎo)向矢量為

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

1.2.4 疼痛評(píng)分 采用11點(diǎn)數(shù)字評(píng)分法，以無痛為“0”，最劇烈疼痛為“10”，共11個(gè)點(diǎn)描述疼痛程度，分值越高，疼痛程度越重?；颊吒鶕?jù)自己的疼痛感覺在相應(yīng)的分值處劃對(duì)號(hào)。

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

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

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

式(23)的共軛復(fù)根的虛部大小可以作為是否是信源的判斷標(biāo)準(zhǔn)，虛部值越接近于0，表示該共軛復(fù)根越有可能是信源對(duì)應(yīng)的共軛復(fù)根，而實(shí)部就是估計(jì)的角度。

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

1)離線計(jì)算獨(dú)立于數(shù)據(jù)的降維矩陣{W1，W2，…，W K};

2)根據(jù)噪聲子空間計(jì)算矢量z;

3)對(duì)第k個(gè)角度區(qū)間，根據(jù)式(22)計(jì)算L+1個(gè)多項(xiàng)式系數(shù)γk1，γk2，…，γk(L+1);

4)對(duì)第k個(gè)角度區(qū)間的多項(xiàng)式式(23)進(jìn)行多項(xiàng)式求根;

5)重復(fù)步驟3和4直到所有角度區(qū)間處理完畢，得到所有的共軛復(fù)根，然后根據(jù)復(fù)根的虛部大小判斷信源對(duì)應(yīng)的根，也就是挑選虛部值最小的P個(gè)根，根據(jù)根的實(shí)部直接計(jì)算信源的角度。

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

3 仿真實(shí)驗(yàn)

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

圖4 本方法獲得的MUSIC譜對(duì)應(yīng)的共軛復(fù)根分布

圖5 均方根誤差隨信噪比的變化關(guān)系

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

圖6 均方根誤差隨快拍數(shù)的變化關(guān)系

4 結(jié)束語(yǔ)

針對(duì)非均勻線陣DOA估計(jì)帶來的高計(jì)算量問題，本文提出了一種快速求根MUSIC方法。該方法利用小角度范圍內(nèi)的導(dǎo)向矢量可以利用低階多項(xiàng)式很好近似的特性，將非均勻線陣DOA估計(jì)問題轉(zhuǎn)化為多組低階實(shí)多項(xiàng)式求根問題。該方法能夠在性能基本保持不變的情況下大大降低標(biāo)準(zhǔn)求根MUSIC算法的計(jì)算復(fù)雜度，非常適合于工程應(yīng)用。

[1]王永良，陳輝，彭應(yīng)寧.空間譜估計(jì)理論與算法[M].北京:清華大學(xué)出版社，2004.

[2]劉周，李軍，劉紅明，等.MIMO雷達(dá)非嚴(yán)格正交信號(hào)單脈沖測(cè)角方法[J].雷達(dá)科學(xué)與技術(shù)，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]林文鳳，周曉青，甘露，等.一種可自動(dòng)配對(duì)的二維波達(dá)方向估計(jì)算法[J].雷達(dá)科學(xué)與技術(shù)，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]黃秋欽，余嘉，梁炎夏，等.基于改進(jìn)PASTd的MUSIC算法的DSP實(shí)現(xiàn)[J].雷達(dá)科學(xué)與技術(shù)，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]康兆敏.數(shù)值計(jì)算方法[M].北京:清華大學(xué)出版社，2013.

[14]蘇峰，王昌海，徐征.基于最小二乘的時(shí)差定位算法[J].雷達(dá)科學(xué)與技術(shù)，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)