基于循環(huán)特性的均勻圓陣快速空間譜計算*

黃秀瓊1，郝克鋼，樊 榮1，萬 群**

(1.中國西南電子技術(shù)研究所，成都610036；2.電子科技大學(xué) 電子工程學(xué)院，成都611731)

基于循環(huán)特性的均勻圓陣快速空間譜計算*

黃秀瓊1，郝克鋼2，樊 榮1，萬 群**2

(1.中國西南電子技術(shù)研究所，成都610036；2.電子科技大學(xué) 電子工程學(xué)院，成都611731)

針對在陣元數(shù)多、分辨率要求高的情況下，基于均勻圓陣的空間譜直接計算方法存在計算量大、實時性差的缺點，利用循環(huán)矩陣、離散傅里葉變換、卷積運算之間的內(nèi)在聯(lián)系，提出了一種快速空間譜計算方法。仿真結(jié)果表明，在保持測向性能完全一致的情況下，所提快速算法降低空間譜計算過程中角度搜索階段的計算量隨陣元數(shù)增加而增加；且在相同計算量條件下，所提快速算方法可以用于更高精度和分辨率的角度搜索。

均勻圓陣；空間譜估計；循環(huán)矩陣；離散傅里葉變換；卷積運算

1 引 言

利用天線陣列計算空間譜是計算信號來波方向、進行波束形成等重要應(yīng)用的關(guān)鍵依據(jù)。在電子偵察、雷達、通信、聲吶等諸多領(lǐng)域的實際應(yīng)用中，要求空間譜計算技術(shù)具有很高的實時性。因此，研究快速的空間譜計算方法具有重要意義。

基于子空間的空間譜計算方法具有超分辨的優(yōu)點，典型的MUSIC算法的計算量主要集中在子空間分解和角度搜索兩部分。過去的幾十年，學(xué)者們在快速子空間分解上開展了大量研究：Xu和Kailath[1]提出的快速子空間分解算法使用Lanczos快速方法得到協(xié)方差矩陣的上三角分解矩陣，然后對其特征分解來獲得信號子空間；Melissions和Tufts[2]提出了由協(xié)方差矩陣的多項式近似信號子空間的方法；Karhunen[3]提出使用離散傅立葉變換(Discrete Fourier Transform,DFT)和余弦變換來近似信號子空間的方法；張柯等人[4]提出一種基于多級維納濾波器(Multi-stage Wiener Filter,MSWF)的矢量陣快速來波方向估計算法，選取矢量陣參考陣元聲壓通道的輸出作為期望信號，通過MSWF的遞推運算得到信號子空間，無需計算陣列協(xié)方差矩陣及特征值分解運算。以上研究工作極大降低了子空間分解的計算量。

在角度搜索階段，需要按照一定的角度間隔在一定的角度范圍內(nèi)確定每一個搜索的角度對應(yīng)的空間譜值，所有空間譜值構(gòu)成空間譜。因此，需要計算的空間譜值的個數(shù)等于要搜索的角度個數(shù)。為了獲得高精度、高分辨率的空間譜計算結(jié)果，需要設(shè)置的角度間隔很小，相應(yīng)需要搜索的角度個數(shù)很多，因而導(dǎo)致計算空間譜所需的時間增加；當(dāng)陣元數(shù)變大時，由于每個搜索的角度所需的計算量變大，又導(dǎo)致計算空間譜所需的時間隨著陣元數(shù)的增大而巨增[5]?？梢?，角度搜索的計算量不容忽視，降低空間譜計算角度搜索階段的計算量，可以提高高精度、高分辨率的空間譜計算方法的計算效率。

目前，針對降低角度搜索階段計算量的研究比較少，而本文利用均勻圓陣導(dǎo)向矢量間循環(huán)移位相等特性與卷積運算的內(nèi)在聯(lián)系有效降低了空間譜計算角度搜索階段的計算量。對算法計算量的分析表明，所提算法降低的計算量與陣元數(shù)呈正比。

2 信號模型

考慮一個具有M個陣元、半徑為R的均勻圓陣，每個陣元所在的角度為

θm=2π(m-1)/M,m=1,2,…,M。

(1)

假設(shè)有來自K個方向的遠場窄帶信號沖激該陣列，則陣列的接收信號為

(2)

式中：s(n)為陣列參考點圓心處接收到的基帶信號，v(n)為接收機噪聲。陣列方向向量表示為

(3)

式中：徑波比α=R/λ，方向向量的第m個元素為am(θk)=ej2παcos(θk-θm)。

本文所研究的方法在空間譜直接計算方法基礎(chǔ)上，利用均勻圓陣導(dǎo)向矢量間循環(huán)移位相等特性與卷積運算的內(nèi)在聯(lián)系，以及離散傅里葉變換與卷積運算之間的對偶關(guān)系，降低空間譜直接計算方法角度搜索階段的計算量，提高空間譜計算效率。

3 空間譜快速計算方法

在噪聲為高斯白噪聲的假設(shè)下，均勻圓陣接收信號的協(xié)方差矩陣譜分解表示為

(4)

式中：λk(k=1,2,…,K)為K個大特征值，其所對應(yīng)的特征向量uk(k=1,2,…,K)構(gòu)成信號子空間的一組標(biāo)準(zhǔn)正交基；λn(n=K+1,K+2,…,M)為M-K個小特征值，其所對應(yīng)的特征向量un(n=K+1,K+2,…,M)構(gòu)成噪聲子空間的一組標(biāo)準(zhǔn)正交基。

本文在已經(jīng)完成子空間分解得到信號子空間基向量的前提下，研究均勻圓陣空間譜角度快速搜索算法。因此，在算法相關(guān)的討論中，不涉及子空間分解?？焖俚淖涌臻g分解算法可以參見文獻[1-4]。

基于子空間的空間譜直接計算方法，對應(yīng)每個搜索角度，通過下式計算其空間譜值[5]：

(5)

將空間譜直接計算方法中的P個搜索角度劃分為M個子集，M為陣元個數(shù)，對應(yīng)地得到M個均勻圓陣方向向量子集為

{a(θq+(m-1)2π/M)}m=1,2,…,M。

(6)

式中：θq=(q-1)·2π/P，q=1,2,…,Q，Q=[P/M]對應(yīng)m=1方向向量子集，[]表示四舍五入取整。

結(jié)合式(3)可以看到，相鄰兩個方向向量子集中對應(yīng)θq的方向向量的元素間存在循環(huán)移位關(guān)系[3]：

(7)

式中：移位矩陣為M階方陣。

用式(6)中各個方向向量子集中對應(yīng)θq的方向向量構(gòu)成M×M濾波矩陣：

(8)

(9)

Wk(θq)的第m個元素為

(10)

上式等號右邊的轉(zhuǎn)置即為式(5)中求和項：

(11)

式中：θp=θq+(m-1)·2π/M。將式(10)、(11)代入式(5)得[6]

(12)

Wk,θq(i)=x(i)?uk(i) 。

(13)

式中：符號?表示循環(huán)卷積，*表示取共軛。

根據(jù)循環(huán)卷積與傅里葉變換的關(guān)系，由關(guān)系式(13)可得

FFT(Wk(θq))=b(θq)×vk。

(14)

式中：符號×表示向量對應(yīng)元素相乘;

符號FFT(·)表示快速傅里葉變換，符號Reverse(·)表示將向量倒序排列。

由(14)式可得

Wk(θq)=IFFT(b(θq)×vk) 。

(15)

在子空間分解得到信號子空間基向量已經(jīng)完成的前提下，現(xiàn)將空間譜快速計算算法整理如下：

Step1 確定角度搜索范圍、相鄰搜索角度之間的間隔d(單位：度)、所有搜索的角度的個數(shù)P。

Step2 將P個搜索角度劃分為M個子集{a(θq+(m-1)2π/M)}，m=1,2,…,M，將第一個子集{a(θq),q=1,2,…,Q}中Q個方向向量倒序排列后進行快速傅里葉變換得{b(θq),q=1,2,…,Q}。

Step3 基于子空間方法求解樣本自相關(guān)矩陣信號子空間的所有列向量{uk,k=1,2,…,K}并對其進行快速傅里葉變換后再取共軛得到{vk,k=1,2,…,K}。

Step4 將集合{b(θq)}中的每個向量分別與集合{vk}的所有的向量元素對應(yīng)相乘再求快速逆傅里葉變換得Q個譜向量集合{Wk(θq),k=1,2,…,K}(q=1,2,…,Q)，每個集合中有K個譜向量。

Step5 將Step 4中得到的Q個集合代入式(12)來計算P個搜索角度上的偽空間譜。

4 算法計算量分析及對比

計算量統(tǒng)計是以一次復(fù)數(shù)乘法和一次復(fù)數(shù)加法為一次運算，即O(1)。設(shè)來波信號個數(shù)為K，陣元個數(shù)為M，搜索方向個數(shù)為P。

空間譜直接計算方法的角度搜索階段，由式(5)可知，每計算一個搜索角度上的空間譜值需要求得K個信號子空間基向量與搜索角度對應(yīng)的陣列方向向量內(nèi)積后再進行取模平方求和，計算量約為O((M+1)K)。總共有P個搜索角度，所以空間譜直接計算方法的總計算量為

O(P(M+1)K) 。

(16)

本文所提空間譜快速計算方法的Step 2中，Q個方向向量的快速傅里葉變換計算量為Q×O(MlbM)；Step 3中，K個信號子空間向量{uk,k=1,2,…,K}的快速傅里葉變換計算量為K×O(MlbM)；Step 4中，基于快速逆傅里葉變換得到Q個譜向量集合的計算量約為Q×K×O(MlbM+M)；Step 5中，式(12)求得所有P個搜索角度上的空間譜值需要進行Q×M×K次復(fù)數(shù)取模平方求和，計算量約為O(QMK)，故總計算量為所有步驟計算量之和：

O(PK(lbM+2)+PlbM+KMlbM) 。

(17)

從式(16)和式(17)可知，空間譜直接計算方法的計算量與陣元數(shù)M呈線性關(guān)系，而本文所提的快速計算方法的計算量與陣元數(shù)的對數(shù)lbM相關(guān)，由此可以推斷出快速算法減少的計算量隨陣元數(shù)增加而增加。換言之，陣元數(shù)越多，快速算法所需的時間越短。為了更為直觀地展示快速算法的有效性，下面給出了兩種算法計算量隨陣元數(shù)和搜索角度數(shù)變化的對比圖。

考慮均勻圓陣獲得高分辨率空間譜的要求，在式(16)和(17)中，設(shè)置搜索角度個數(shù)P=720；考慮角度估計個數(shù)上限為陣元數(shù)M，信號來波方向個數(shù)設(shè)為K=?M/2」(符號?·」表示向下取整)，約為可估計角度個數(shù)上限的一半。圖1給出了空間譜直接計算方法和空間譜快速計算方法計算量隨陣元個數(shù)M的變化。

圖1 角度搜索計算量隨陣元數(shù)變化(P=720,K=?M/2」 )Fig.1 Computational complexity in angle searching phase vs. the number of elements(P=720,K=?M/2」 )

從圖 1中可以看出，在分辨率為0.5°、來波信號方向個數(shù)約為陣元數(shù)一半的條件下，空間譜快速計算方法減少的計算量隨著陣元數(shù)的增加而增加，且陣元數(shù)越多，減少的計算量越多。陣元數(shù)M=12時，本文所提算法將角度搜索階段的計算量減少了52%。

圖2給出了K=6、M=12時，空間譜直接計算方法和空間譜快速計算方法計算量隨搜索角度間隔d=360/P的變化。

圖2 角度搜索計算量隨分辨力變化(M=12,K=6)Fig.2 Computational complexity in angle searching phase vs. the resolution(M=12,K=6)

從圖2中可以看出，相同計算量條件下，空間譜快速搜索算法可以用于更高精度和分辨率的角度搜索。

5 性能一致性驗證

由第3部分的理論推導(dǎo)和式(12)可知，空間譜直接計算方法和本文提出的快速計算方法計算的空間譜完全相同。為了進一步說明所提快速算法在減少計算量的同時保持了空間譜計算的性能，下面給出將兩種算法計算得到的空間譜用于DOA估計的性能對比，如圖3所示。仿真中，參數(shù)設(shè)置與第4節(jié)一致，均勻圓陣陣元數(shù)M=12，信號來波方向個數(shù)K=6，搜索角度個數(shù)P=720。為了得到合理的陣列結(jié)構(gòu)，保證DOA估計性能，設(shè)置陣列徑波比α=4。接收信號測量快拍數(shù)L=20。

圖3 角度估計性能曲線圖(M=12,P=720，K=6)Fig.3 Performance of the DOA estimation(M=12,P=720，K=6)

從圖3中可以看出，兩種方法的空間譜估計性能曲線重合，說明本文提出的快速方法降低計算量是以不損失估計性能為前提的。

6 結(jié)論

在保持測向性能完全一致的情況下，本文所提的空間譜快速計算方法通過降低空間譜計算過程中角度搜索階段的計算量來提高算法的整體計算效率。在相同的計算成本約束下，本文所提的空間譜快速計算方法可用于更高精度和分辨率的角度搜索。

[1] XU G,KAILATH T. Fast subspace decomposition[J]. IEEE Transactions on Signal Processing,1994,42(3):539-551.

[2] TUFTS D,MELISSINOS C D. Simple, effective computation of principal eigenvectors and their eigenvalues and application to high-resolution estimation of frequencies[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1986, 34(10): 1046-1053.

[3] KARHUNEN J T,JOUTSENALO J.Sinusoidal frequency estimation by signal subspace approximation[J].IEEE Transactions on Acoustics, Speech, and Signal Processing, 1992, 40(12): 2961-2972.

[4] 張柯，程菊明，付進. 基于多級維納濾波器的聲矢量陣空間譜估計算法[J].兵工學(xué)報，2015，36(11)：2128-2134.

ZHANG Ke,CHENG Juming,FU Jin.DOA estimation of acoustic vector sensor array based on multi-stage Wiener filter[J].Acta Armamentarii,2015,36(11):2128-2134.(in Chinese)

[5] 齊崇英.適用于均勻圓陣的陣列互耦校正與信源測向方法:101149429[P].2008-03-26.

[6] 萬群.采用均勻圓陣快速測定空間譜的方法: 201410100251.7[P].2014-03-18.[7] BABAK. An acceleration of FFT-based algorithms for the match-count problem[J].Information Processing Letters,2017,125:1-4.

[8] ABEYSEKERAS S. Computationally efficient DOA estimation using sparse arrays with I/Q mismatch and D.C. offsets[C]//Proceedings of 2016 IEEE International Conference on Digital Signal Processing (DSP).Beijng:IEEE,2016:49-53.

[9] BELLILI F, AMOR S B, AFFES S, et al. Low-complexity DOA estimation from short data snapshots for ULA systems using the annihilating filter technique[J/OL].EURASIP Journal on Advances in Signal Processing, 2017,2017(1):1-16[2017-09-20]. https://link.springer.com/content/pdf/10.1186%2Fs13634-017-0480-1.pdf. https://doi.org/10.1186/s13634-017-0480-1.

[10] WEI L, SHAO W, QI W D, et al. Peak-to-peak search: fast and accurate DOA estimation method for arbitrary non-uniform linear array[J].Electronics Letters,2015,51(25):2078-2080.

FastSpatialSpectrumEstimationBasedonCycleCharacteristicofUniformCircularArray

HUANG Xiuqiong1，HAO Kegang2，F(xiàn)AN Rong1，WAN Qun2

(1.Southwest China Institute of Electronic Technology,Chengdu 610036,China； 2.School of Electronic Engineering,University of Electronic Science and Technology of China,Chengdu 611731,China)

In the case of the large number of array elements and high resolution requirement, the algorithm for computing the spatial spectrum directly has low eficiency and poor realtime performance.For this disadvantage,this paper proposes an efficient uniform circular array (UCA) based spatial spectrum estimation algorithm to ruduce the computational complexity in angle searching phase by using the inner relationship between the eigen-decomposition of the cyclic matrix, the discrete Fourier transform (DFT) and the convolution operation.When the performance is kept exactly the same as that of the direct algorithm, the proposed algorithm is able to decrease more computational complexity if there are more array elements. The results of numerical simulation show that the proposed fast algorithm can achieve higher resolution with the same computing resources.

uniform circular array(UCA);spatial spectrum estimation;cyclic matrix;discrete Fourier transform(DFT);convolution operation

10.3969/j.issn.1001-893x.2017.12.009

黃秀瓊，郝克鋼，樊榮，等.基于循環(huán)特性的均勻圓陣快速空間譜計算[J].電訊技術(shù)，2017,57(12):1399-1403.[HUANG Xiuqiong，HAO Kegang，F(xiàn)AN Rong，et al.Fast spatial spectrum estimation based on cycle characteristic of uniform circular array[J].Telecommunication Engineering,2017,57(12):1399-1403.]

2017-10-03；

2017-12-08

date:2017-10-03；Revised date：2017-12-08

國家科技重大專項(2016ZX03001022)；國家自然科學(xué)基金資助項目(U1533125)

wanqun@uestc.edu.cnCorrespondingauthorwanqun@uestc.edu.cn

TN911.7

A

1001-893X(2017)12-1399-05

黃秀瓊(1971—)，女，重慶萬州人，高級工程師，主要從事機載綜合射頻傳感器系統(tǒng)總體技術(shù)研究；

Email: xiuqhuang@sina.com

郝克鋼(1993—)，男，四川雅安人，博士研究生，主要研究方向為空時自適應(yīng)處理技術(shù)；

樊榮(1984—)，男，四川閬中人，博士，工程師，主要從事綜合射頻認知系統(tǒng)設(shè)計與開發(fā)；

Email: fanrong@alu.uestc.edu.cn

萬群(1971—)，男，江西南昌人，教授、博士生導(dǎo)師，主要研究方向為陣列信號處理、認知無線電定位等。

Email: wanqun@uestc.edu.cn