黃秀瓊1，郝克鋼，樊 榮1，萬(wàn) 群**

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

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

黃秀瓊1，郝克鋼2，樊 榮1，萬(wàn) 群**2

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

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

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

1 引 言

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

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

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

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

2 信號(hào)模型

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

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

(1)

假設(shè)有來(lái)自K個(gè)方向的遠(yuǎn)場(chǎng)窄帶信號(hào)沖激該陣列，則陣列的接收信號(hào)為

(2)

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

(3)

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

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

3 空間譜快速計(jì)算方法

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

(4)

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

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

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

(5)

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

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

(6)

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

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

(7)

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

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

(8)

(9)

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

(10)

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

(11)

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

(12)

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

(13)

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

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

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

(14)

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

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

由(14)式可得

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

(15)

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

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

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

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

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

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

4 算法計(jì)算量分析及對(duì)比

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

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

O(P(M+1)K) 。

(16)

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

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

(17)

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

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

圖1 角度搜索計(jì)算量隨陣元數(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°、來(lái)波信號(hào)方向個(gè)數(shù)約為陣元數(shù)一半的條件下，空間譜快速計(jì)算方法減少的計(jì)算量隨著陣元數(shù)的增加而增加，且陣元數(shù)越多，減少的計(jì)算量越多。陣元數(shù)M=12時(shí)，本文所提算法將角度搜索階段的計(jì)算量減少了52%。

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

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

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

5 性能一致性驗(yàn)證

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

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

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

6 結(jié)論

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

[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ìn). 基于多級(jí)維納濾波器的聲矢量陣空間譜估計(jì)算法[J].兵工學(xué)報(bào)，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] 齊崇英.適用于均勻圓陣的陣列互耦校正與信源測(cè)向方法:101149429[P].2008-03-26.

[6] 萬(wàn)群.采用均勻圓陣快速測(cè)定空間譜的方法: 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ì)算[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

國(guó)家科技重大專(zhuān)項(xiàng)(2016ZX03001022)；國(guó)家自然科學(xué)基金資助項(xiàng)目(U1533125)

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

TN911.7

A

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

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

Email: xiuqhuang@sina.com

郝克鋼(1993—)，男，四川雅安人，博士研究生，主要研究方向?yàn)榭諘r(shí)自適應(yīng)處理技術(shù)；

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

Email: fanrong@alu.uestc.edu.cn

萬(wàn)群(1971—)，男，江西南昌人，教授、博士生導(dǎo)師，主要研究方向?yàn)殛嚵行盘?hào)處理、認(rèn)知無(wú)線(xiàn)電定位等。

Email: wanqun@uestc.edu.cn