基于平行因子四線性分解的二維角度和頻率聯合估計

許凌云 張小飛 許宗澤

(南京航空航天大學電子信息工程學院 南京 210016)

1 引言

信號多維參數聯合估計是近年來陣列信號處理發展的一個重要方向，特別是空間到達角和信號源頻率的聯合估計成為一個研究的熱點。大多數的時空參數聯合估計算法都是在1維算法的基礎上擴展而得，如2維MUSIC算法[1,2]、各種ESPRIT算法及2維傳播算子(PM)算法等， MUSIC方法要求譜峰搜索，其多維搜索過程的計算量很大，工程上不可行，而且當信噪比或快拍數下降到一定門限時，其性能迅速下降。文獻[3]利用TLS-ESPRIT算法，結合均勻矩形陣列對2維到達角和頻率進行聯合估計，主要優點在于計算上避免了譜峰搜索，對互譜矩陣進行特征值分解或對接收信號進行奇異值分解來估計參數，但對于面陣/雙平行線陣來說，特征值分解或奇異值分解的運算量很大，并存在參數配對問題，在低噪比情況下參數配對往往出現模糊現象。文獻[4]提出基于傳播算子的2維角度估計方法，與ESPRIT算法相比，計算量大大減小，且估計的參數能夠自動配對，但存在孔徑損失，影響估計精度。

平行因子技術是高維數據分析的一種方法，平行因子概念開始于心理學領域，后來發展到化學、光譜學、色譜學等領域，近年來在信號處理和通信領域，平行因子技術迅速發展，已經逐步成為通信信號處理領域里一種新的研究手段[5?9]。以上文獻討論的都是基于平行因子三線性分解的技術，相比而言，國內外在平行因子4維數據分解辨識領域研究較少。本文給出雙平行線陣列結構，提出了一種基于平行因子四線性分解的2維角度和頻率聯合估計算法。根據陣列輸出數據構建一個四線性模型，通過對此模型分解可以得到2維角度和頻率估計的參數。這種方法是迭代算法，不需要特征值分解或奇異值分解，不需要參數配對，只需一定的迭代次數就能收斂。

文中用到的符號和算子說明如下：diag(·)是對角化算子； ⊙ 表示 Khatri-Rao積； [·]?表示矩陣的Moore-Renrose逆矩陣；||||F表示F范數；Di(·)表示提取矩陣第i行元素構成的對角陣；黑體小寫字母表示矢量；黑體大寫字母表示矩陣。

2 數據模型

圖1 雙平行線陣結構

雙平行線陣的結構如圖1所示，它由兩個均勻平行線子陣1和子陣2組成，每個線陣有M個陣元，以原點為參考點，X軸方向陣元間距為Δx，Y軸方向陣元間距為Δy，設 Δx=Δy=d,d為入射信號最小波長的一半，假設遠場有K個彼此獨立且到達角均不同的窄帶信號入射到該系統上，第k個信號的2維到達角和頻率分別為(θk, φk,fk),k=1,2,… ,K，其中θk,φk和fk分別代表第k個信源的仰角、方位角和頻率。各陣元輸出的噪聲為統計獨立、均值為零和方差為 σ2的加性高斯白噪聲，且噪聲與信號相互獨立。則兩個子陣接收的信號分別為

式中

整個陣列接收的信號為

為了對頻率進行估計，對陣列接收信號加上時延 τp(p=1,2,… ,P)后采樣，可得到

假設空時信道在L個符號內恒定不變，定義可表示為

定義

將2 × (P+ 1)個切片累積成M× (P+1)×2×L的4維數據集，則式(6)有以下四線性模型形式[10]：

其中am,k是Ax矩陣中第(m,k)元素，bn,k是Ay矩陣中第(n,k)元素，ψp,k是矩陣Φ中第(p,k)個元素，sl,k是S中第(l,k)個元素，wn,m,p,l形成了4維噪聲數據集。式(7)對xn,m,p,l的分解稱為 4維平行因子分析，也稱為四線性分解。根據其對稱性，可寫成其它3種切片形式：

⊙表示Khatri-Rao積，同理式(8)，式(9)和式(10)也可以表示成2ML×(P+ 1)的矩陣Y,2L(P+ 1)×M的矩陣Z和ML(P+ 1)× 2的矩陣U，如式(12)-式(14)所示。

從上可知雙平行線陣中2維角度和頻率的估計可轉化為對方向矩陣Ax,Ay和時延矩陣Φ3個矩陣的估計。

3 基于平行因子四線性分解的聯合估計算法

3.1 平行因子四線性模型分解

四線性交替最小二乘(QALS)算法是四線性模型數據檢測的一種常用方法。本文采用QALS算法對式(11)給出的平行因子模型進行擬合。迭代過程中的代價函數為

給定Ax,Ay和Φ,S的最小二乘解為

利用式(12)，給定S,Ax和Ay,Φ的最小二乘解為

利用式(13)，給定Φ,S和Ay,Ax的最小二乘解為

利用式(14)，給定Φ,S和Ax,Ay的最小二乘解為

利用上述方法，每次更新一個矩陣，本次更新的矩陣立刻參與下一個矩陣的求解，直至算法收斂。

3.2 可辨識性分析

四線性模型的本質特征就是其唯一性，首先說明k秩的概念，它在多線性代數中起著非常重要的作用。

定義 1 對于給定的矩陣A∈?I×F，當且僅當A包含至少r個但不包含r+1個線性獨立的列時，A的秩為rA=rank(A)=r。如果矩陣A的任意k列獨立，則A的k秩kA=k，矩陣k秩比矩陣秩的要求更為嚴格，且kA≤rA。當矩陣A列滿秩時，kA=rA。

定理1[10]=1,2;p=1,2,… ,P+ 1)， 其 中A∈?M×K,A∈xy?2×K,Φ∈?(P+1)×K,S∈?K×L，如果矩陣Ax,Ay,Φ,S的k秩滿足下面的條件

則Ax,Ay,Φ和S對于列交換和(復數)尺度變換是唯一的，即Ax,Ay,Φ,S可唯一確定。

3.3 頻率和2維角度估計的實施步驟

上述算法結束后，得到時延矩陣Φ的估計值，設矩陣的第k列為。在沒有噪聲的情況下，對其取相角可以得到

angle(·)表示取元素的相角運算，其中τ=[0,2πτ1,…, 2 πτ]T，令P，其中1P+1表示P+1個全1矢量，則

fk的最小二乘結果為

根據式(21)-式(23)，可依次得到K個信源的頻率估計。

m的最小二乘結果為

同理可得到，即的估計值。所以2維DOA估計為

算法總的步驟可總結如下：

步驟 1 經過四線性最小二乘算法迭代運算后得到Ax,Ay和Φ的估計值和，并令和分別為3個估計矩陣的第k列矢量；

步驟2 按照式(21)-式(23)可依次得到K個信源的頻率估計；

步驟3 按照式(24)-式(29)可依次得到K個信源的2維角度估計。

本文算法的復雜度主要集中在四線性最小二乘的迭代運算上，每次迭代計算復雜度(復乘)為O(4K3+8KLM(P+1))及O(K2[2M+ML+L(P+1)+ 2(P+1) + 2M(P+1) + 2LM+ML(P+1) + 2L(P+1)])，文獻[8]采用的三線性分解的算法運算量主要在三線性最小二乘的迭代運算上，每次計算復雜度為O(3K3+6KM(P+1)L+K2[2M(P+1)+L(P+1)+2ML])，迭代的次數n跟參數估計的數目和初始值的設置有關。文獻[4]提出的ESPRIT算法集中在協方差的計算、三次特征值的分解及配對搜索上，算法運算量為O{(2PMK)2L+(2PMK)3+4(P?1)MK2+2MK2+ 2(M?1)K2+ 6K3}。

4 計算機仿真

為了驗證文中提出的算法的性能，采用 Monte Carlo仿真來評估算法的2維角度和頻率估計性能。假設噪聲為加性高斯白噪聲，仿真中采用圖1所示的雙平行線陣列。假設非相關窄帶信號源數為K=3，其2維DOA和頻率分別為(15o,10o, 1.2 MHz), (25o,20o, 1.5 MHz), (35o,30o, 1.8 MHz)，陣元間距 Δx=Δy=75 m，快拍數L=500，定義RMSE：RMSE=，其中amk為第k個信源的第m個Monte Carlo仿真的角度或頻率的估計值，aok為第k個信源的角度或頻率的真實值。

仿真 1 仿真采用陣元數為8×2的雙平行線陣列，進行了50次Monte Carlo仿真，數據長度為L=500。圖2給出了本文算法在SNR＝10 dB時方位角和仰角的散布圖及方位角和頻率的散布圖，從圖2可以看出此算法具有較好的聯合角度和頻率估計性能。

圖2 SNR=10 dB的2維角度和頻率散布圖

仿真 2 為了說明本文算法在不同信噪比下角度和頻率參數聯合估計性能，仿真參數同上，Monte Carlo試驗次數為 500，數據長度為L=500。圖3-圖11中橫坐標 SNR的定義為 10lg(PS/PN)，其中PS為信號平均功率，PN為噪聲平均功率，為了進行比較，ESPRIT算法和基于三線性(trilinear)分解算法的結果一起在圖3-圖5中給出。圖3-圖5分別為在不同信噪比條件下的頻率和2維DOA估計值的RMSE曲線，可以看出，隨著信噪比的提高，RMSE變小，估計性能更好。同時可看出本文算法的性能優于其它兩種算法，trilinear分解又要優于ESPRIT算法。這主要因為三線性和四線性分解的算法充分利用接收端信息的同時，其每次迭代都有最小的二乘閉式解，而四線性分解參與最小二乘擬合的數據增多，提高了參數估計的精度。文獻[10]也從理論上進一步證明了四線性分解的 CRB要小于三線性分解的CRB。而ESPRIT算法在特征值分解和配對搜索過程中存在誤差積累，參數估計精度下降。

圖3 仰角隨信噪比變化性能曲線

圖4 方位角隨信噪比變化性能曲線

圖5 頻率隨信噪比變化性能曲線

仿真 3 本文為了說明在不同快拍數下頻率和2維角度聯合估計的性能，采用3個非相關窄帶信號源，快拍數分別設為10, 50, 100, 300, 500。圖6-圖8分別給出了角度和頻率在不同快拍數下隨信噪比變化的RMSE曲線，從圖中可以看出，在快拍數很小的情況下，性能較好，隨著快拍數的增加，參數估計性能越來越好。從仿真中可知，該算法在小樣本數下也能較好的工作。

圖6 仰角在不同快拍數下性能曲線

圖7 方位角在不同快拍數下性能曲線

圖8 頻率在不同快拍數下性能曲線

仿真 4 本文為了說明在不同天線數下頻率和2維角度聯合估計的性能，采用3個非相關窄帶信號源， 陣元數分別設為6×2,8×2,10×2,12×2。圖9-圖11分別給出了角度和頻率在不同天線數下隨信噪比變化的RMSE曲線，從圖中可以看出，隨著天線數的增加，參數估計性能更好，表明該算法有較好的空間分集性能。

圖9 仰角在不同天線數下性能曲線

圖10 方位角在不同天線數下性能曲線

圖11 頻率在不同天線數下性能曲線

5 結束語

本文將平行因子四線性分解技術應用到雙平行線陣列中提出一種新的2維角度和頻率估計算法。該算法首先利用四線性交替最小二乘算法估計出方向矩陣和頻率矩陣，然后利用頻率矩陣的Vandermonde特征和方向矩陣本身的結構特點及最小二乘方法進行運算可得到頻率和2維角度的估計值。該方法無需譜峰搜索和參數配對，與基于三線性分解的算法和 ESPRIT算法進行了比較和分析，具有更高的估計精度，并且在小樣本數下也能較好地工作，更利于實際應用。仿真結果驗證了該方法的有效性。

[1] Swindlehurst A L, Stoica P, and Jansson M. Exploiting arrays with multiple invariances using music and mode[J].IEEE Transactions on Signal Processing, 2001, 49(11):2511-2521.

[2] Wax M, Shan T J, and Kailath T. Spatio-temporal spectral analysis by eigenstructure methods[J].IEEE Transactions on Acoustics,Speech,Signal Processing, 1984, 32(4): 817-827.

[3] Strobach P. Total least squares phased averaging and 3-D ESPRIT for joint azimuth-elevation carrier estimation[J].IEEE Transactions on Signal Processing, 2001, 49(1): 54-62.

[4] 任勛立, 廖桂生, 曾操. 一種低復雜度的二維波達方向估計方法[J]. 電波科學學報, 2005, 20(4): 526-530.Ren Xun-li, Liao Gui-sheng, and Zeng Cao. A low complexity method for estimating 2D-DOA[J].Chinese Journal of Radio Science, 2005, 20(4): 526-530.

[5] Sidiropoulos N D, Bro R, and Giannakis G B. Parallel factor analysis in sensor array processing[J].IEEE Transactions on Signal Processing, 2000, 48(8): 2377-2388.

[6] 張劍云, 鄭志東, 李小波. 雙基地雷達收發角及多普勒頻率的聯合估計算法[J]. 電子與信息學報, 2010, 32(8): 1843-1848.Zhang Jian-yun, Zheng Zhi-dong, and Li Xiao-bo. An algorithm for DOD-DOA and Doppler frequency jointly estimating of bistatic MIMO radar[J].Journal of Electronics&Information Technology, 2010, 32(8): 1843-1848.

[7] Zhang X, Gao X, and Xu D. Novel blind carrier frequency offset estimation for OFDM system with multiple antennas[J].IEEE Transactions on Wireless Communications, 2010, 9(3):881-885.

[8] Zhang Xiao-fei, Feng Gao-peng, and Xu Da-zhuang. Anglefrequency estimation using trilinear decomposition of the oversampled output[J].Wireless Personal Communications,2009, 51(2): 365-373.

[9] 劉旭, 許宗澤, 雷磊. 應用矩陣分解的陣列信號參數辨識[J].應用科學學報, 2010, 28(1): 49-55.Liu Xu, Xu Zong-ze, and Lei Lei. Identification of array signal parameters based on matrix decomposition[J].Journal of Applied Science, 2010, 28(1): 49-55.

[10] Liu Xiang-qian and Sidirodpoulos N D. Cramer-Rao lower bounds for low-rank decomposition of multidimensional arrays[J].IEEE Transactions on Signal Processing, 2001,49(9): 2074-2086.