求根MUSIC初值設置和更新算法

閆鋒剛，張 薇，金 銘

(哈爾濱工業(yè)大學(威海)信息與電氣工程學院，264209山東威海)

信號波達方向(direction-of-arrival，DOA)估計［1-2］是雷達、聲吶、無線通信和無源定位等應用中經(jīng)常遇到的重要研究課題.以多重信號分類(multiple signal classification，MUSIC)［3］和旋轉不變子空間(estimate signal parameters via rotational invariance technique，ESPRIT)［4］為代表的子空間類算法的提出，實現(xiàn)了傳統(tǒng)空間譜估計向超分辨測角的飛躍，但 MUSIC算法龐大的計算量和ESPRIT算法較低的估計精度阻礙了超分辨算法的工程化進度［5］.

在均勻線陣(uniform linear array，ULA)下，陣列導向矢量具有特殊的范德蒙結構，使得MUSIC算法繁雜的譜峰搜索簡化為多項式求根問題，由此促成了求根 MUSIC(root-MUSIC)算法［5］的誕生.雖然 root-MUSIC僅是 MUSIC在ULA下的特例，然而理論分析和實驗結果證明:root-MUSIC的估計精度比MUSIC的估計精度更高［6-7］，因而其相比于MUSIC實現(xiàn)了低復雜度和高精度的雙贏.文獻［8］通過理論分析論述了root-MUSIC估計性能與環(huán)境、陣元參數(shù)的內(nèi)在聯(lián)系;文獻［9］借助中心對稱陣列結構實現(xiàn)了root-MUSIC的實值化，在降低計算量的同時提高了算法的估計精度;文獻［10-11］借助流型分離(manifold separation technique，MST)技術，將任意陣列的導向矢量近似地分裂成兩個獨立部分的乘積，利用包含信號角度信息且具有范德蒙結構的部分、以多項式求根實現(xiàn)信號 DOA估計.文獻［12］利用任意陣列結構下的MUSIC空間譜是關于信號DOA周期函數(shù)的事實，提出了傅里葉域求根MUSIC(fourier-domain，F(xiàn)D-MUSIC)算法，將root-MUSIC由ULA推廣到了任意陣列結構.文獻［13］基于傅里葉級數(shù)的勞倫結構，對 FD-root-MUSIC進一步進行了改進，有效降低了多項式的求根階數(shù)，降低了算法的復雜度.

上述算法豐富了root-MUSIC的理論內(nèi)涵，但它們在實際工程應用中一般均需借助迭代算法對多項式求根.選擇合適的初值并基于一定的原則對其更新和控制以使迭代快速、準確地收斂具有重要研究意義.

本文在分析root-MUSIC多項式根分布特點的基礎上，提出了一種迭代初值到單位圓平均距離最短(least average distance to unit circle，LADTUC)的準則，進而基于該準則發(fā)展了一種適用于root-MUSIC的初值設置和更新算法.理論分析和試驗結果表明，該算法能有效避免迭代過程中的錯誤解并加速迭代收斂速度，從而為root-MUSIC的實際工程化提供理論參考.

1 數(shù)據(jù)模型及root-MUSIC算法原理

1.1 數(shù)據(jù)模型

考察xoy平面上由M個以半波長間隔放置的陣元組成的ULA，假設陣列各通道附加高斯白噪聲且空間有L個輻射源，定義波達方向DOA為信號來向與陣列法線的夾角，則陣列一次快拍接收數(shù)據(jù)為［3］

式中:x(t)是M×1維接收數(shù)據(jù)向量;s(t)是L×1

維信號向量;n(t)是M×1維噪聲向量.A(θ)是M×L維導向矢量矩陣且導向矢量a(θ)定義為

陣列接收數(shù)據(jù)協(xié)方差矩陣定義為

式中:RS為信號協(xié)方差矩陣，σ2n為噪聲功率.對R進行特征值分解，可得

式中:ρi，i=1，2，…，L和ρj，j=L+1，…，M分別為R的L個大特征值及(M-L)個小特征值，ei和ej分別為是ρi和ρj對應的特征矢量.定義:

則S和G的列分別張成信號子空間和噪聲子空間.

1.2 root-MUSIC算法原理

依據(jù)上述陣列和數(shù)據(jù)模型，root-MUSIC算法定義如下所示的多項式［5］:

其中p(z)定義為

對f(z)求根即可獲得DOA.由于f(z)存在z*項，求根較為復雜.為便于計算，做如下修正:

此時，多項式f(z)的階數(shù)為2(M-1)，所以有(M-1)對根，每一對根是相互共軛的關系.這些根中的L個根z1，z2，…，zL正好分布在單位圓上，且

實際中，得到多項式f(z)中L個接近于單位圓的根后，可按下式求解信號DOA:

2 LADTUC準則及最佳初值設置

2.1 f(z)最小值設置

通常為得到DOA，只需根據(jù)式(8)令f(z)=0并求根即可，但當快拍數(shù)或信噪比較低時，在信號入射方向有

因此，為得到精確DOA，需估計fmin并對修正多項式f(z)-fmin=0進行求根.

由MUSIC零譜統(tǒng)計性能知，發(fā)射能量為Ei(i=1，2，…，L)的輻射源，在 MUSIC 空間零譜上產(chǎn)生的最小值為［14-15］

式中J為快拍數(shù)，輻射源能量Ei(i=1，2，…，L)在一般情況下是難以獲得的，一種近似的求解方法是認為每個輻射源的輻射能量都相同，即

式中:E為觀測信號總能量，ρj是協(xié)方差矩陣R特征值分解后對應的噪聲空間特征值，其數(shù)值等于噪聲的能量.由上述分析，可設置

實際中，當快拍數(shù)或信噪比較高時，fmin接近于0.此時，可直接選取近似值fmin≈0.

2.2 最佳迭代初值設置

對方程f(z)-fmin=0采用牛頓法求解［16］，得

式中:f'(z)為f(z)的一階導函數(shù)，z0為迭代初值，zn是經(jīng)過第n次迭代后得到的解.

可選取一定容差?，當z0與zn滿足

時迭代結束，然后將zn帶入式(11)即可求得DOA.

迭代初值z0的選取至關重要，其直接決定著迭代的收斂速度和正確性.由于代表信號DOA的根分布在單位圓上，故z0的選擇應以到單位圓平均距離最短為準則，這便是本文提出的LADTUC準則，后續(xù)仿真實驗證實了該準則的合理性.

由于DOA取值范圍為［-π/2，π/2］，只需考慮半單位圓即可.為了敘述便利，以下均用U代表上半單位圓弧.設R(cosθ，sinθ)為U上任意一點，其中，θ為R與圓心連線的傾角，顯然，θ符合均勻分布特性，其概率密度函數(shù)為

設最佳迭代初值為Q1，則Q1到R(cosθ，sinθ)平均距離的平方L21為

依據(jù)LADTUC原則，令

易求得，最佳迭代初值為

2.3 迭代初值更新方法

LADTUC原則保證的是平均意義下的最佳，實際中的某次迭代解算可能因受到U內(nèi)、外根的干擾而無法收斂到U上.為了解決該問題，定義數(shù)列Qn:

式中:點A為弧Un，n=1，2，…上任意一點，Un為將U進行k1倍等分后、順時針方向的第k2+1份弧.易求得k1、k2與n滿足如下關系:

式中符號「·?表示向上取整.

由式(21)可見:Qn，n=1，2，… 是 LADTUC意義下對Q1的擴展，其代表的是到弧Un平均距離最小的點列.設弧Un上、下端點和圓心連線的傾角為θ1和θ2，易求得

利用類似于式(18)、(19)的方法，可得Qn，n=1，2，… 的橫、縱坐標值為

顯然，當n→ ∞ 時，θ1→θ2→0，Qn→1.因此，Qn的所有元素均在U內(nèi).圖1給出了Qn的前7個元素在U內(nèi)的相對位置關系.基于Qn，可對迭代初值更新，以保證迭代最終收斂于U上.

圖1 數(shù)列Q的前7個元素

假設U上分布正確根R1和R2，U內(nèi)、外分別分布干擾根Ⅰ1、Ⅰ2和O1、O2，如圖 2 所示.依據(jù)LADTUC原則，首先設置迭代初值z0為Q1，若迭代收斂于干擾根Ⅰ1，則更新z0為Ⅰ1相對于Q1反方向的Q2，以避免Ⅰ1的干擾;否則，更新z0為Q3.若更新后迭代仍收斂于干擾根Ⅰ2，則更新z0為Ⅰ2相對于Q3反方向的Q5，以避免Ⅰ2的干擾;同理，當?shù)俅问諗坑诟蓴_根O1時，則更新z0為O1相對于Q5反方向的Q4，以避免O1的干擾，….如此反復更新z0，直至迭代收斂于正確根R1.

在LADTUC原則下，每次初值更新都是以平均最接近U的等分弧段而進行.為防止多次更新仍不收斂，可預設最大容忍迭代次數(shù)N，當?shù)螖?shù)達到N/2時，下一步則應“逃逸”至相反的半弧進行.如圖2所示，當設置z0為Q4仍不收斂時，下一步應更新z0為Q3，以在U的右半弧尋求正確根R2.

圖2 LADTUC原則下迭代初值設置和更新

3 收斂性能分析

式(8)根的總數(shù)為

其中，代表輻射源DOA的根有L對，因此，干擾根的個數(shù)為

采用LATDUC原則進行迭代初值設置和更新時，當干擾根距離數(shù)列Q第k個元素Qk的距離小于Qk到Uk的平均距離時，則迭代會收斂到干擾根而不會收斂到單位圓上，此時收斂會失敗.因而，第k次初值設置收斂失敗的概率為

因此，經(jīng)過k次初值更新，成功收斂概率為

由于0＜qk＜1，由式(28)可見:pLADTUC隨k的增加急劇增大，因而本文算法每經(jīng)一次初值更新后，成功收斂的概率都會大幅增加.

若采用隨機法設置迭代初值，只進行1次初值設置迭代收斂成功(收斂到U上)的概率為

由于隨機法的每次初值設置相互獨立，故其經(jīng)過k次初值設置成功收斂的概率恒為prand，這意味著隨機設置和更新初值的收斂概率并不會隨k的增加而增大.因而，LADTUC相比于隨機設置、更新迭代初值法從統(tǒng)計平均上提高了成功收斂的概率.

4 仿真試驗

試驗采用ULA陣列，陣元間距為半個波長.設快拍數(shù)J=256，迭代停止容差?=10-2，最大容忍迭代次數(shù)N=30.設輻射能量彼此相等，信源DOA通過隨機生成［-π/2，π/2］上均勻分布的L個角度得到.Monte Carlo試驗次數(shù)均為500次.

4.1 最佳迭代初值Q1合理性驗證

本試驗用于驗證式(20)給出最佳迭代初值Q1(0，2/π)的合理性.試驗設置陣元數(shù)M=10，迭代初值為z0=0+Δi，其中步長以間隔0.05從0.05增大到1.試驗分為如下兩個部分.

首先，選取L=3并以RSN為參變量，統(tǒng)計Δ的取值與采用牛頓法求解方程根的迭代次數(shù)及收斂概率的關系，結果如圖3所示.試驗中，當?shù)偞螖?shù)超過N=30時認為收斂失敗.

圖3 信噪比參變時，步長與迭代次數(shù)及收斂概率關系

其次，選取RSN=5 dB并以信號數(shù)L為參變量，統(tǒng)計步長值與采用牛頓法求解方程根的迭代次數(shù)及收斂概率的關系，結果見圖4.試驗中，當?shù)偞螖?shù)超過N=30時則認為收斂失敗.

由圖3、4可見:相比于RSN，L對迭代次數(shù)和迭代的成功收斂概率影響更強，這是由于RSN僅影響方程(8)根的分布，而L除了影響方程根分布之外還決定著干擾根的多少.另外，由圖3、4可觀察到:當Δ=0.65≈2/π時，成功收斂概率最大，迭代次數(shù)也接近最小，這證明了初值Q1(0，2/π)的合理性.

4.2 LADTUC初值更新原則的合理性驗證

本試驗用于驗證圖2所示的基于LADTUC原則進行迭代初值更新方法的合理性.試驗設置陣元數(shù)M=15，RSN=5 dB，然后分別按LADUC原則和隨機法更新迭代初值，并按如下兩方面進行平均結果統(tǒng)計.

首先，統(tǒng)計迭代初值更新次數(shù)與信號數(shù)L的關系，結果見圖5.試驗中，當?shù)偞螖?shù)超過N=30時則認為收斂失敗.

圖4 信號數(shù)參變時，步長與迭代次數(shù)及收斂概率關系

圖5 迭代初值更新次數(shù)與信號數(shù)目關系

其次，設置初值更新的最大容忍次數(shù)為3，然后統(tǒng)計迭代成功收斂概率與信號數(shù)L的關系，結果見圖6.試驗中，當?shù)偞螖?shù)超過N=30或初值更新次數(shù)超過3次時，均認為收斂失敗.

由圖5可見:基于LADTUC原則進行迭代初值更新，最多只需要進行3次以下的更新即可成功收斂，其相比于隨機更新法顯著減少了更新次數(shù).另一方面，當信號數(shù)L較小時，本文方法所需的更新次數(shù)比隨機法的更新次數(shù)更少.這是因為L越小時方程(8)的干擾根較多，而LADTUC原則相比于隨機法從平均意義上保證了迭代收斂概率更高，因此其所需的初值更新次數(shù)越少.

由圖6可見，當L＞2時，LADTUC原則下的迭代成功收斂概率均＞95%;而隨機法的成功收斂概率僅為85%左右，其相比于本文方法低10個百分點.與圖5相似，從圖6可見:本文方法相比于隨機發(fā)在信號數(shù)L較少(即式(8)干擾根較多)時的優(yōu)勢更為明顯，迭代收斂概率更高.上述實驗結果說明了依據(jù)LADTUC準則對root-MUSIC算法進行初值設置和更新是合理性的.

圖6 迭代成功收斂概率與信號數(shù)目關系

5 結語

本文在分析root-MUSIC算法根分布特點的基礎上，提出了一種利用工程實現(xiàn)的root-MUSIC初值設置和更新算法.該算法充分利用了LADTUC原則對于迭代收斂的控制作用，保證了root-MUSIC算法迭代收斂的快速性和正確性，為root-MUSIC算法的實際工程化提供了一定的理論參考.

［1］JIN H.Joint space-time parameter estimation for multicarrier CDMA systems ［J］.IEEE Transactions on Vehicular Technology，2012，61(7):3306-3311.

［2］KHABBAZIBASMENJ A.Efficient transmit beamspace design for search-free based DOA estimation in MIMO radar［J］.IEEE Transactions on Signal Processing，2014，62，(6):1490-1500.

［3］YAN F G，JIN M，LIU S，et al.Real-valued MUSIC for efficient direction estimation with arbitrary array geometries ［J］. IEEE Transactions on Signal Processing，2014，62(6):1548-1560.

［4］ZHENG G M，CAO B X，YANG M L.Unitary ESPRIT algorithm for bistatic MIMO radar［J］.Electronics Letters，2012，48(3):1122-1123.

［5］CHENG Q，HUANG L，So H C.Improved unitary root-MUSIC for DOA estimation based on pseudo-noise resampling［J］.IEEE Signal Processing Letters，2014，21(2):140-144.

［6］STOCIA P，NEHORAI A.MUSIC，maximum likelihood，and Cramer-Rao bound:further results and comparisons［J］.IEEE Transactions on Acoustics，Speech and Signal Processing，1990，38(12):2140-2150.

［7］CAO Y，ZHANG Z，DAI F，et al.Direction of arrival estimation for monostatic multiple-input multiple-output radar with arbitrary array structures ［J］.IET Radar，Sonar＆ Navigation，2012，6(7):679-686.

［8］LI F，LIU H，VACCARO R J.Performance analysis for DOA estimation algorithms:Unification，Simplification，and Observations［J］.IEEE Transactions on Aerospace and Electronic Systems，1993，29(4):1170-1184.

［9］MARIUS P，ALEX B G，MARTIN H.Unitary root-MUSIC with a real-valued eigendecomposition:a theoretical and experimental performance study［J］.IEEE Transactions on Signal Processing，2000，48(5):1306-1314.

［10］BELLONI F，RICHTER A，KOIVUNEN V.Extension of root-MUSIC to non-ULA array configurations［C］//International Conference on Acoustics，Speech and Signal Processing(ICASSP).Toulouse，F(xiàn)rance:IEEE，2006:4(5):14-19.

［11］BELLONI F，RICHTER A，KOIVUNEN V.DOA estimation via manifold separation for arbitrary array structures ［J］. IEEE Transactions on Signal Processing，2007，55(10):4800-4810.

［12］RUBSAMEN M，GERSHMAN A B.Direction-of-arrival estimation for non-uniform sensor arrays:from manifold separation to Fourier domain MUSIC methods［J］.IEEE Transactions on Signal Processing，2009，57(2):588-599.

［13］ZHUANG J，LI W，MANIKAS A.Fast root-MUSIC for arbitrary arrays ［J］.Electronics Letters，2010，46(2):174-176.

［14］HARRY B，MICHAEL L，WENGROVITZ S.Statistical characterization of the MUSIC null spectrum［J］.IEEE Transactions on Signal Processing，1991，39(6):1911-1921.

［15］BARABELL A J.Improving the resolution performance ofeigenstructure based direction finding algorithms［C］//International Conference on Acoustics，Speech and Signal Processing(ICASSP).Boston，America:IEEE，1983:336-339.

［16］CADZOW J A.Multiple source location-the signal subspace approach ［J］. IEEE Transactions on Acoustics， Speech and SignalProcessing， 1990，38(7):1110-1125.