基于數(shù)據(jù)插值的二維MUSIC譜峰搜索算法

馬英英　劉帥　金銘

摘 要： 針對(duì)MUSIC算法在進(jìn)行二維到達(dá)角估計(jì)時(shí)運(yùn)算量大、估計(jì)速度較慢的問題，對(duì)基于數(shù)據(jù)插值二維MUSIC譜峰搜索算法進(jìn)行研究。在MUSIC算法的基礎(chǔ)上先用較大步長在二維空間進(jìn)行粗略搜索，之后在粗搜譜峰位置附近的小區(qū)間內(nèi)進(jìn)行二維插值，最后用小步長在該區(qū)間內(nèi)進(jìn)行譜峰搜索，得到精確的角度估計(jì)值。仿真實(shí)驗(yàn)表明，該算法可有效地減少二維空間搜索的運(yùn)算量。

關(guān)鍵詞： MUSIC算法； 譜峰搜索； 運(yùn)算量； 數(shù)據(jù)插值

中圖分類號(hào)： TN911.7?34 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 1004?373X（2015）15?0035?04

Two?dimensional MUSIC spectrum?peak searching algorithm based on data interpolation

MA Yingying， LIU Shuai， JIN Ming

（School of Information and Electrical Engineering， Harbin Institute of Technology （Weihai Branch）， Weihai 264209， China）

Abstract： Since multiple signals classification （MUSIC） algorithm has large calculation quantity and low estimation speed when estimating two?dimensional angle of arrival （AOA）， two?dimensional MUSIC spectrum?peak searching algorithm based on data interpolation is studied. On the basis of MUSIC algorithm， larger step is applied to proceed rough searching in two?dimensional space， then two?dimensional interpolation is conducted in minizone near searched spectrum?peak position. Spectrum?peak searching is carried on in the zone with smaller step， and precise angle estimation value is obtained. Simulation experiments indicate that the algorithm can reduce calculation quantity in two?dimensional space searching effectively.

Keywords： MUSIC algorithm； spectrum?peak searching； calculation quantity； data interpolation

0 引 言

空間譜估計(jì)測向技術(shù)[1]作為一門新興的空間信號(hào)處理技術(shù)，在雷達(dá)、聲納、氣象、通信等相關(guān)領(lǐng)域中都有著廣泛應(yīng)用[2] 。多重信號(hào)分類[3]（Multiple Signal Classification，MUSIC）算法作為一種高分辨率的譜估計(jì)方法使得空間譜估計(jì)技術(shù)進(jìn)入了超分辨測角階段。經(jīng)典MUSIC算法因?yàn)樾枰M(jìn)行協(xié)方差矩陣的計(jì)算、特征值分解和譜峰搜索，計(jì)算量很大，特別是在二維或多維空間進(jìn)行搜索時(shí)，計(jì)算量尤其大，很難在工程中應(yīng)用。后來人們在其基礎(chǔ)上提出了很多改進(jìn)算法。其中，1985年Roy和Kailath提出的借助旋轉(zhuǎn)不變技術(shù)的參數(shù)估計(jì)算法[4]（Estimating Signal Parameters Via Rotational Invariance Techniques，ESPRIT）是利用陣列流形的某些特性形成一個(gè)可以直接求解的函數(shù)，能夠比較方便地得到所需要的估計(jì)參數(shù)；Barabell提出的求根MUSIC算法[5]根據(jù)多項(xiàng)式求根進(jìn)行參數(shù)估計(jì)。但上述兩種算法對(duì)陣列形式均有要求，求根MUSIC算法僅能用于一維陣列且ESPRIT算法在二維條件下存在參數(shù)配對(duì)問題。而另一部分改進(jìn)算法還是不能避免譜峰搜索[6?8]，如基于空間平滑技術(shù)去相關(guān)的MUSIC（MMUSIC）算法[9]，在譜峰搜索上的運(yùn)算量還是很大。因此有學(xué)者提出，可以結(jié)合數(shù)值方法中的插值法進(jìn)行譜峰搜索，可有效地減少運(yùn)算量。

數(shù)據(jù)插值方法[10]是數(shù)值分析中的最基本方法之一，主要解決的問題是根據(jù)離散數(shù)據(jù)構(gòu)造一個(gè)簡單易于計(jì)算的函數(shù)代替原有的復(fù)雜函數(shù)。本文將插值算法應(yīng)用在譜峰搜索當(dāng)中，可以大大減少搜索的點(diǎn)數(shù)，從而減少計(jì)算量，并且對(duì)陣列形式無特殊要求。

1 MUSIC算法二維空間譜估計(jì)

用如圖1所示的均勻平面陣來估計(jì)二維空間譜的方位角與俯仰角。

設(shè)陣元位置為[（xk，yk）， k=1，2，…，M，]信號(hào)入射參數(shù)為[（θi，φi）， i=1，2，…，N，]分別表示方位角與俯仰角，其中方位角表示與[x]軸的夾角。

以原點(diǎn)處為參考陣元，則各陣元相對(duì)參考陣元的時(shí)延為：

[τki=1c（xkcosθicosφi+yksinθicosφi）] （1）

圖1 均勻平面陣

對(duì)于遠(yuǎn)場窄帶信號(hào)，整個(gè)天線陣列所接收到的信號(hào)為：

[X（t）=AS（t）+N（t）] （2）

式中：[X（t）]為陣列的[M×1]維快拍數(shù)據(jù)矢量；[N（t）]為陣列的[M×1]維噪聲數(shù)據(jù)矢量；[S（t）]為空間信號(hào)的[N×1]維矢量；[A]為空間陣列的[M×N]維流型導(dǎo)向矢量陣，且：

[A=[a1（φ，θ） a2（φ，θ） … aN（φ，θ）]] （3）

其中，導(dǎo)向矢量：

[ai（φ，θ）=ej2πλ（x1sinθicosφi+y1sinθisinφi），…，ej2πλ（xNsinθicosφi+yNsinθisinφi）T] （4）

式中[λ]為波長。

得到陣列的接收信號(hào)模型后進(jìn)行二維空間譜估計(jì)。首先求陣列輸出的相關(guān)矩陣[R∈CN×N]，有：

[R=1Li=1LXXH] （5）

對(duì)[R]進(jìn)行特征分解，則其[N-K]個(gè)較小特征值對(duì)應(yīng)的特征向量張成了[R]的噪聲子空間[UN]。利用噪聲子空間與信號(hào)方向矢量正交即[UN⊥span{A}]的特性，得空間二維MUSIC算法定義的空間譜函數(shù)為：

[PMUSIC=1aH（φ，θ）UNUHNa（φ，θ）] （6）

對(duì)空間譜函數(shù)進(jìn)行二維搜索，獲得的譜峰處即指示信號(hào)來波方向的[（φk，θk）]。

2 插值算法

數(shù)學(xué)上常用的插值法[11]有Lagrange插值，Newton插值，Hermite插值等多項(xiàng)式插值法以及分段低次插值法，分段三次樣條插值等分段插值法。多項(xiàng)式插值法在插值多項(xiàng)式次數(shù)較高時(shí)會(huì)出現(xiàn)不收斂的龍格現(xiàn)象；分段低次多項(xiàng)式插值通常計(jì)算簡單但曲線光滑性差；分段三次多項(xiàng)式插值要求插值函數(shù)具有連續(xù)的二階導(dǎo)數(shù)，可以得到具有更高光滑度的曲線。鑒于MUSIC譜圖的譜峰非常尖銳，所以為得到更高的精度本文選用分段三次樣條插值法，并與分段線性插值法，分段三次多項(xiàng)式插值法作比較。

2.1 分段線性插值

對(duì)給定區(qū)間[[a，b]]做劃分[Δ：a=x0

（1） [φ（x）]在區(qū)間[[a，b]]上連續(xù)；

（2） [φ（xi）=yi，i=0，1，2，…，n]；

（3） [φ（x）]在每個(gè)子區(qū)間[[xi，xi+1]]上是線性函數(shù)。

則稱[φ（x）]為[[a，b]]上關(guān)于數(shù)據(jù)[（xi，yi），i=0，1，2，…，n]的分段線性插值函數(shù)。

同理，若[φ（x）]在每個(gè)子區(qū)間[[xi，xi+1]]上是三次多項(xiàng)式，則稱[φ（x）]為[[a，b]]上關(guān)于數(shù)據(jù)[（xi，yi），i=0，1，2，…，n]的分段三次多項(xiàng)式插值函數(shù)。

2.2 分段三次樣條插值

對(duì)給定區(qū)間[[a，b]]做劃分[Δ：a=x0

（1）[S（x）]在[[a，b]]上具有連續(xù)的二階導(dǎo)數(shù)；

（2）[S（xi）=yi，i=0，1，2，…，n]；

（3）[S（x）]在每個(gè)子段[[xi，xi+1]]上都是三次多項(xiàng)式。

則稱[S（x）]為分段三次樣條函數(shù)。

當(dāng)[x∈[xi，xi+1]]時(shí)，記[hi=xi+1-xi，][S（x）]的一般表達(dá)式為：

[S（x）=Mi（xi+1-x）6h3+Mi+1（x-xi）36h+yi-Mih26 xi+1-xh+yi+1-Mi+1h26 x-xih] （7）

式中：[Mi，i=0，1，2，…，n]是需要待定的系數(shù)，且[Mi=][S（xi）]。

據(jù)[S（xi+0）=S（xi-0）]，可得[n-1]個(gè)條件：

[μiMi-1+2Mi+λiMi+1=di，i=1，2，…，n-1] （8）

[式中：][λi=hi+1hi+hi+1，][μi=hihi+hi+1=1-λi，][di=6hi+hi+1?] [yi+1-yihi+1-yi-yi-1hi=][6f（xi-1，xi，xi+1）]。。。

由于[n+1]個(gè)數(shù)據(jù)點(diǎn)需要[n+1]個(gè)系數(shù)，而式（8）只提供了[n-1]個(gè)條件，所以通過端點(diǎn)約束，即增加邊界條件來補(bǔ)充另外兩個(gè)條件。邊界條件一共有三種，分別是：

（1） 第一種邊界條件

[S（x0）=f（x0），S（xn）=f（xn）] （9）

（2） 第二種邊界條件

[S（x0）=f（x0），S（xn）=f（xn）] （10）

若[S（x0）=S（xn）=0，]則稱之為自然邊界，滿足自然邊界的樣條函數(shù)稱為自然樣條，它是通過所有數(shù)據(jù)點(diǎn)的插值函數(shù)中總曲率最小的惟一函數(shù)，因此自然三次樣條是插值所有數(shù)據(jù)點(diǎn)的最光滑函數(shù)。

（3） 第三種邊界條件

當(dāng)被插值函數(shù)是以[b-a]為周期的周期函數(shù)時(shí)，則要求[S（x）]也是周期函數(shù)，這時(shí)，邊界條件為：

[S（x0+0）=S（xn-0），S（x0+0）=S（xn-0），S（x0+0）=S（xn-0）] （11）

這樣就可以解出[n+1]個(gè)系數(shù)了，再將其代入到式（7）就可以得到相應(yīng)子區(qū)間內(nèi)的樣條函數(shù)。

在利用上述插值方法對(duì)MUSIC算法粗搜的二維空間譜進(jìn)行插值時(shí)，先利用得到的數(shù)據(jù)在[x]軸方向進(jìn)行插值，再用同樣的插值法對(duì)[y]軸方向進(jìn)行插值，這樣就得到了經(jīng)二維插值后的MUSIC空間譜圖。

3 仿真及分析

實(shí)驗(yàn)中采用8陣元[L]陣，信號(hào)入射角為（30.45°，60.27°），在高斯白噪聲條件下進(jìn)行。信噪比為30 dB，快拍數(shù)為1 024。性能分析實(shí)驗(yàn)中蒙特卡洛次數(shù)為100次。仿真結(jié)果的插值區(qū)域均為以粗搜譜峰為中心，邊長為2°的正方形區(qū)域。

3.1 基于插值的譜峰搜索算法實(shí)現(xiàn)

圖2為搜索步長為0.1°時(shí)，MUSIC算法在二維空間中搜索的譜圖；圖3為在MUSIC算法0.1°搜索的基礎(chǔ)上，用分段線性插值法進(jìn)行插值后的譜圖；圖4為在MUSIC算法0.1°搜索的基礎(chǔ)上，用分段三次多項(xiàng)式插值法進(jìn)行插值后的譜圖；圖5為在MUSIC算法0.1°搜索的基礎(chǔ)上，用分段三次樣條插值法進(jìn)行插值后的譜圖?？梢钥闯鼋?jīng)插值后的譜圖要比圖2密集得多，理論上譜峰搜索的精度也會(huì)提高，下面通過仿真來比較幾種算法的性能。

圖2 二維MUSIC算法0.1°粗搜譜圖

圖3 分段線性插值譜圖

3.2 性能比較

圖6，圖7為基于插值算法的二維MUSIC譜峰搜索算法均方根誤差（RMSE）比較分析。實(shí)驗(yàn)中將MUSIC算法0.1°搜索及MUSIC算法0.01°搜索的RMSE與基于分段線性插值、分段三次多項(xiàng)式插值及分段三次樣條插值的二維MUSIC譜峰搜索算法的RMSE進(jìn)行比較，分別得到方位角與俯仰角的RMSE對(duì)比曲線。結(jié)果表明：分段線性插值搜索的RMSE相比于MUSIC算法0.1°粗搜并沒有改善，即分段線性插值在二維譜峰搜索中并不適用；分段三次多項(xiàng)式插值搜索與分段三次樣條插值搜索的RMSE相比于MUSIC算法0.1°粗搜都有改善；分段三次樣條插值搜索的性能要比分段三次多項(xiàng)式插值搜索更好且與MUSIC算法0.01°搜索的RMSE相差很小。

圖4 分段三次插值譜圖

圖5 分段三次樣條插值譜圖

圖6 方位角RMSE隨信噪比的變化

3.3 運(yùn)算速度分析

表1為用兩種算法進(jìn)行譜峰搜索實(shí)驗(yàn)100次的時(shí)間，可以看到基于三次樣條插值的二維MUSIC譜峰搜索的方法明顯減少了運(yùn)算量，提高了二維空間譜峰搜索的速度。

圖7 俯仰角RMSE隨信噪比的變化

表1 算法耗時(shí)比較

[算法＼&MUSIC算法0.01°搜索＼&分段三次樣條插值＼&耗時(shí) /s＼&151.755＼&41.536＼&]

4 結(jié) 論

MUSIC算法作為一種高分辨率的譜估計(jì)方法使得空間譜估計(jì)技術(shù)進(jìn)入了超分辨測角階段，但譜峰搜索的過程運(yùn)算量很大，尤其是在二維空間中。所以如何快速準(zhǔn)確地進(jìn)行譜峰搜索是當(dāng)前需要解決的重要問題之一，本文主要研究了基于多種插值算法的MUSIC譜峰搜索算法。

在二維空間中建立了平面[L]陣的陣列信號(hào)模型，仿真得到了二維插值搜索的空間譜圖，通過對(duì)算法性能在參數(shù)影響和運(yùn)算速度兩個(gè)方面的對(duì)比分析，得到如下結(jié)論：

（1） 基于分段線性插值的二維MUSIC譜峰搜索算法的RMSE相比于MUSIC算法0.1°粗搜并沒有提高，所以分段線性插值法在二維MUSIC譜峰搜索中不適用。

（2） 基于分段三次多項(xiàng)式插值及分段三次樣條插值的二維MUSIC譜峰搜索的RMSE相比于MUSIC算法0.1°粗搜都有改善；分段三次樣條插值搜索的性能要比分段三次多項(xiàng)式插值搜索更好且與MUSIC算法0.01°搜索的RMSE相差很小。

（3） 基于分段三次樣條插值的二維MUSIC譜峰搜索算法相比MUSIC算法0.01°搜索明顯減少了運(yùn)算量，大大提高了二維空間譜峰搜索的速度。

參考文獻(xiàn)

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

[2] KRIM H， VIBERG M. Two decades of array signal processing research [J]. IEEE Signal Processing Magazine， 1996， 13（4）： 67?69.

[3] SCHMIDT R O. Multiple emitter location and signal parameter estimation [J]. IEEE Transactions on Antennas and Propagation， 1986， 34（3）： 276?280.

[4] GAO F， GERSHMAN A B. A generalized ESPRIT approach to direction?of?arrival estimation [J]. Signal Processing Letters， 2005， 12（3）： 254?257.

[5] RAOB D， HARI K V S. Performance analysis of root?music. IEEE Transactions Acoustics， Speech and Signal Processing， 1989， 37（12）： 1939?1949.

[6] 吳世龍，羅景青，劉有軍.一種用于二維MUSIC算法的譜峰識(shí)別方法研究[J].信號(hào)處理，2005，21（z1）：16?19.

[7] 吳江華，周圍.DOA估計(jì)的一種改進(jìn)MUSIC算法[J].無線電通信技術(shù)，2008，34（1）：39?42.

[8] 周小軍，譚薇，馮大政，等.基于解相干的改進(jìn)MUSIC算法DOA估計(jì)[J].無線電工程，2014，44（12）：18?21.

[9] 何子述，黃振興.修正MUSIC算法對(duì)相關(guān)信號(hào)源的DOA估計(jì)性能[J].通信學(xué)報(bào)，2000，21（10）：14?17.

[10] MATHEWS J H，F(xiàn)INK K D.數(shù)值方法（Matlab版）[M].北京：電子工業(yè)出版社，2010.

[11] 張麗娟.三種插值方法的應(yīng)用與比較[J].赤峰學(xué)院學(xué)報(bào)：自然科學(xué)版，2010，26（3）：1?3.