基于張量分解的子空間頻率估計及應用

摘 要:為了提高短時序列的頻率估計精度，提出一種基于張量分解的適用于多通道接收數據的子空間高分辨率的頻率估計算法。該方法將接收到的多通道采樣信號按照一定的形式構建三維數據結構，利用子空間的移不變性，采用張量的高階奇異值分解直接對數據進行頻率估計。仿真表明在低信噪比的情況下，基于張量分解的頻率估計算法性能要優于基于矩陣奇異值分解的頻率估計算法，將此算法應用于目標軌跡測量中，可以獲得高精度的目標軌跡參數。關鍵詞:頻率估計; 張量分解; 奇異值分解;多普勒頻率

中圖分類號:TP873文獻標識碼:A

文章編號:1004-373X(2010)18-0138-03

Subspace Frequency Estimation and Application Based on Tensor Decomposition

HAN Feng， WEI Guo-hua， WU Si-liang

(School of Information and Electronics， Beijing Institute of Technology， Beijing 100081， China)

Abstract: Subspace high-resolution algorithm of frequency estimation suitable for multi-channel receiving data is proposed to improve the accuracy of short-sequence frequency estimation based on the tensor decomposition. The sampling signal received by multi-channel are structured into three-dimensional data structure according to certain form. The frequency may be estimated by using the character of subspace invariance and high-order singular value decomposition. The simulation results show that the performance of the frequency estimation algorithm based on tensor decomposition is better than the one based on matrix singular value decomposition algorithm in the case of the low SNR. The algorithm may obtain high accurate target trajectory parameters， which is applied to target trajectory measurement.Keywords: frequency estimation; tensor decomposition; singular value decomposition; Doppler frequency

0 引 言

在短時序列的頻率估計中，經典的譜估計算法存在瑞利限的限制，估計精度低;而基于子空間的高分辨率譜估計算法如ESPRIT算法，則可以取得較高的頻率估計精度[1-2]，基于子空間分解的頻率估計就是在含寬帶噪聲的觀測中抽取有用信息的過程，數據矩陣的奇異值分解可以用來揭開隱藏的低秩信息，利用奇異值或特征值的大小，可以將信號子空間分為信號子空間和噪聲子空間兩部分，特征子空間方法就是用信號或噪聲子空間進行低秩信息的提取。子空間分解分為2類，一類是基于矩陣的奇異值分解(SVD);另一類是基于張量的高階奇異值分解(HOSVD)。基于矩陣的SVD分解是強有力的工具，已被廣泛應用到信息獲取，協作濾波，計算語言學，計算機視覺等領域[3]。然而SVD分解被限制在二維的數據處理領域，在一些應用中，數據具有三維或更高的模式，如在多傳感器技術、陣列信號處理中，高維數據的獲取越來越容易，與低維數據相比包含了更多的信息，高維數據分析已成為信號與信息學領域的一大熱點，為了分析高維數據、得到更加緊湊的數據模型，準確估計信號參數，需要一種可以對高維數據的進行分解的工具，基于張量的HOSVD分解則可以滿足這種要求[4]。

張量分解包括多種形式，如INDSCAL，PARAFAC2，CANDELING，PARATUCK2等，其中有2類特別的分解形式被認為是矩陣奇異值分解(SVD)的高階推廣:CANDECOMP/PARAFAC(CP)分解是將張量分解為一階張量和的形式;Tucker分解是主元素分析的高階形式(PCA)。Lathauwer等證明了這兩種方法是高維數據的廣義多重線性SVD，由于數據維數可以大于三維，稱為高階奇異值分解。信號處理團體正不斷地嘗試引入HOSVD作為信號處理的工具，一方面使多重線性代數能夠成為譜分析的一個強有力的分析工具，另一方面通過在譜估計問題的分析能夠促進對張量概念的發展，數值計算和子空間方法的研究，使這些算法能夠更好地適用于目標運動參數估計，核磁共振(MRS)，多通道腦電信號(EEG)等實際問題[4-7]。

1 張量分解原理

張量是多重線性分析中的基本對象，也是對向量、矩陣等概念的擴展，向量與矩陣都是張量的特例，即一維的張量叫向量，二維的張量叫矩陣，三維的張量叫盒(box)。在張量中，維(dimension)稱作模式(mode)，把張量A沿n模式展開為矩陣的操作叫展開(unfold)，記為A(n)。Tucker分解是高階主成分分析的一種形式，它將張量分解成一個核函數在每個模上乘以矩陣的形式，如三維張量A∈iI×J×K的分解如下所示:

A≈G×1A(1)×2A(2)×3A(3)=

∑Pp=1∑Qq=1∑Rr=1gpqrap#8226;aq#8226;ar=G;A(1)，A(2)，A(3)(1)



式中:A(1)∈iI×P;A(2)∈iJ×Q;A(3)∈iK×R稱為矩陣因子，通常是正交的，可以認為是每個模中的主成分，G被稱為核張量。高階奇異值分解示意圖如圖1所示。

圖1 高階奇異值分解示意圖

張量A中的每個元素可以表示為:

xijk≈∑Pp=1∑Qq=1∑Rr=1gpqraipbjqckr(2)

式中:i=1，2，…，I;j=1，2，…，J;k=1，2，…，K;P，Q，R稱為矩陣因子的主成分數，分解算法可以理解為:A(i)提供了沿著第i個模式的最大能量方向，由于求A(n)的方法就是把A沿第i個模式展開，進行奇異值分解后的左奇異值矩陣(其列矢量為特征矢量)，因此N階張量的HOSVD就是計算N個I1I2…In…IN不同矩陣的SVD。

三階張量的R階秩估計是對于給定的復數三階張量，找到一個復數三階秩(R1×R2×R3)張量，使目標函數f()最小。

f()=‖A-‖2(3)

在矩陣運算中，最優秩估計方法可以通過對矩陣SVD截斷獲得，張量的最優秩，即，(R1，R2，R3)估計是通過HOSVD的截斷獲取。最優秩的估計方法可以通過子空間的處理計算得到，Tucker-ALS的秩:(R1，R2，R3)估計算法描述如下:

procedure HOSVD(X，R1，R2，R3)

for n=1，…，N do

A(n)←取X(n)的Rn個左奇異向量

end for

G←X×1A(1)×2A(2)×1A(3)

return G，A(1)，A(2)，A(3)

end procedure

關于張量分解的原理、性質和算法已有多篇文獻闡述，本文不再詳細推導[5，8]。

2 算法描述

假設Q為接收通道數目，一次快拍數據包含N×Q個復數采樣點，在沒有噪聲的情況下，多通道的數據模型表示為指數信號和模型:

x(q)n=∑Kk=1Akexp{jψ(q)k}exp{jωkn}=∑Kk=1c(q)kznk(4)

式中:K為信號分量數;n=0，1，2，…，N-1，表示一次快拍中每個通道采集的數據;q為通道號，1≤q≤Q;ck:ck=Akexp{jψk}，第k個信號分量的復數幅度和相位;ψ(q)k為通道的相位;zk為第k個信號分量的極化形式，zk=exp(jωk);ωk為多普勒數字角頻率，ωk=2πfd/fs，fd是信號頻率;fs是數字抽樣頻率。

三維張量的構成方式如圖2所示。每個通道接收到的數據都可以構成一個三維張量H，張量中的元素可表達如下:

hi1i2i3=x(i1-1)+(i2-1)+(i3-1)=xi1+i2+i3-3

圖2 三維張量構成方式

式中:1≤i1≤I1，1≤i2≤I2，1≤i3≤I3，I1+I2+I3=N+2，每通道的數據堆疊在一起，形成三階張量的塊形式。

hi1i2i3=∑Kk=1c(i3)k(zi1-1k#8226;zi2-1k#8226;zi3-1k)(5)

H=∑Kk=1ck#8226;1z1kz2kzI1k#8226;1z1kz2kzI2k#8226;1z1kz2kzI3k

=G×1A(1)×2A(2)×3A(3)(6)

式中:

A(1)=1…1

z11…z1B

z21…z2B



zI11…zI1B，A(2)=1…1

z11…z1B

z21…z2B



zI21…zI2B，

A(3)=1…1

z11…z1B

z21…z2B



zI31…zI3B

對張量分解，得到最優秩:(K，K，K)估計，對矩陣A(n)通過總體最小二乘方法可以得到特征值λk，因為:

λk=k=exp{jωk}=exp{j2πfkfs}

所以:

fk=Im[ln(λk)]#8226;fs/2π

3 頻率的估計精度

設信號包含2個頻率分量，f1=2 000 Hz，f2=2 200 Hz，信號抽樣頻率為80 kHz，幀長度為190，分別采用SVD和HOSVD方法估計信號頻率，仿真次數為500次，頻率的均方根誤差如圖3所示。仿真中可以看出，在信噪比低于10 dB時，基于張量的頻率估計精度要優于基于矩陣SVD分解算法的頻率估計精度。

圖3 SVD算法與HOSVD算法的參數估計精度統計

4 算法在目標運動參數估計中的應用

在雷達陣列信號處理中，常采用短基線天線陣列接收非合作運動目標的回波信號，提取出目標回波中的多普勒頻率用于目標軌跡的提取[9-10]，假設點目標以速度v做勻速直線運動，目標運動軌跡的偏角為α，目標運動軌跡的傾角為β，坐標原點到運動軌跡直線的距離為r0，運動軌跡上，目標初始位置到最近點的距離為L，接收通道數目為Q，目標運動參數的估計問題可以歸結非線性最小優化問題:

0，，^0=argminr0，v，L0∑N0n=1d(tn)-L0-vtnr20+(L0-vtn)22

式中:tn為離散的觀測時間;N0為觀測的快拍次數。

設雷達發射波長為0.1 m，基帶信號采樣率為80 kHz，目標以勻速直線運動從300 m處飛向脫靶點。采用均方根誤差(RMSE)來評價估計參數的性能。回波信號是在理想信號的基礎上加入信噪比為10 dB高斯白噪聲，接收通道Q=8，表1給出經過500次蒙特卡洛仿真，不同運動軌跡參數的估計精度，可以看出參數r0的小于0.3 m，速度的估計精度小于真值的5‰。

表1 目標運動參數的估計結果統計

r0 /mv /(m/s)

第1組真值6300

RMSE0.060.08

第2組真值12320

RMSE0.070.05

第3組

真值18400

RMSE0.040.05

第4組

真值24500

RMSE0.130.07

第5組

真值30600

RMSE0.260.10

基于SVD的方法，計算量與數據段的長度有關，設幀長為N，則矩陣分解需要與N3成正比次數的乘法運算，而HOSVD方法除了與數據段的長度有關外，還與秩估計算法的迭代次數有關，高階張量正交迭代算法(HOOI)，它是交替最小二乘(ALS)的一種，收斂速度呈線性，為了減小誤差距離，需要進行多次迭代;Rayleigh商迭代法(RQI)的收斂速度是平方收斂，迭代次數越多，計算量越大，目前，新的梯度算法的引入，如基于BFGS擬牛頓的方法，可加速最優秩估計算法的收斂，且多核并行計算技術的發展，使得基于HOSVD方法的目標參數估計能夠在目標捕獲后5 min內給出參數估計結果。

5 結 語

基于張量分解的頻率估計算法的原理是子空間的旋轉不變性質和信號模型的物理意義，實質是通過增加噪聲噪聲子空間和信號子空間中奇異值的距離提高參數估計的精度，數據仿真表明這種可用于多通道數據的頻率估計方法可以改善頻率的估計精度。

參考文獻

[1]PETRE S，RANDOLPH M.現代信號譜分析[M].吳仁彪，韓萍，馮青，等譯.北京:電子工業出版社，2007.

[2]MANOLAKIS D G， INGLE V K， KOGNE S M. Statistical and adaptive signal processing[M]. London: Artech House Inc.， 2005.

[3]魏木生.廣義最小二乘問題的理論和計算[M].北京:科學出版社，2006.

[4]TAMARA G K， BRETT W B. Tensor decompositions and applications[J]. SIAM Review， 2009， 51(3): 455-500.

[5]NION D， LATHAUWER L De. A block component model-based blind DS-CDMA receiver[J]. IEEE Transactions on Signal Processing， 2008， 56(11): 5567-5579.

[6]LATHAUWER L De， CASTAING J. Tensor-based techniques for the blind separation of DS-CDMA signal[J]. IEEE Transactions on Signal Processing， 2007， 87(2): 322-336.

[7]HAARDT M， ROEMER F， DEL G G. Higher-order SVD-based subspace estimation to improve the parameter estimation accuracy in multidimensional harmonic retrieval problems[J]. IEEE Transactions on Signal Processing， 2008， 56(7): 3198-3213.

[8]BOYER R， LATHAUWER L De， ABED-MERAIM K. Higher order tensor-based method for delayed exponential fitting[J]. IEEE Transactions on Signal Processing， 2007， 55(6): 2795-2809.

[9]魏國華，吳嗣亮，王菊，等.脫靶量測量技術綜述[J].系統工程與電子技術，2004，26(6):768-772.

[10]馬淑芬.矢量脫靶量測量系統數據處理方法的研究[D].北京:北京理工大學，2007.