劉魯濤 曹 瑩* 鄭 昱

①(哈爾濱工程大學(xué)信息與通信工程學(xué)院 哈爾濱 150001)

②(南京電子技術(shù)研究所 南京 210039)

1 引言

波達(dá)方向(Direction Of Arrival, DOA)估計是戰(zhàn)術(shù)和戰(zhàn)略應(yīng)用中許多電子戰(zhàn)的主要組成部分，近代以來，由于軍事和民用的迫切需求，DOA估計發(fā)展十分迅速[1–4]。多重信號分類(MUltiple SIgnal Classification, MUSIC)算法[5]作為經(jīng)典的DOA估計算法，具有精度高、分辨力高和普遍的陣列適用性等優(yōu)點(diǎn)。將普通陣元用極化敏感陣元代替之后，并按照特定的順序在空間中排列的陣列系統(tǒng)稱為極化敏感陣列[6]。和普通陣列相比，極化敏感陣列具有更高的抗干擾能力、抗角模糊能力、信號檢測能力等優(yōu)點(diǎn)[7]。因此，對基于極化敏感陣列的MUSIC算法的硬件實(shí)現(xiàn)的研究具有重要意義。

目前，DOA估計多采用數(shù)字信號處理器(Digital Signal Processor, DSP)多核運(yùn)算實(shí)現(xiàn)，但DSP為串行操作，使得基于MUSIC算法的DOA估計局限在毫秒級，無法滿足DOA估計的實(shí)時性。隨著FPGA(Field Programmable Gate Array)技術(shù)的迅速發(fā)展，F(xiàn)PGA芯片的容量越來越大，提供的知識產(chǎn)權(quán)(Intellectual Property, IP)庫也越來越全面，F(xiàn)PGA并行度高的優(yōu)勢逐漸展露，為提高DOA估計的實(shí)時性提供了更好的硬件選擇。文獻(xiàn)[8]提出了天線陣列縮減和協(xié)方差歸一化技術(shù)，在100 MHz的工作頻率下僅耗時2 μs即可完成DOA估計，但其提出的技術(shù)僅適用于常規(guī)的標(biāo)量陣，并不適用于極化敏感陣列，且其使用代數(shù)計算的特征分解方法僅適用于協(xié)方差矩陣維數(shù)小于5的情況，有其局限性。文獻(xiàn)[9]實(shí)現(xiàn)MUSIC算法的2維搜索，使用了雅各比(Jacobi)算法的串行脈動陣列結(jié)構(gòu)以及多尺度峰值所搜方法，在100 MHz的工作頻率下完成1次DOA估計需要1 ms。文獻(xiàn)[10]提出了一種完全實(shí)數(shù)計算的MUSIC算法，為了縮減并行Jacobi算法的資源消耗問題，提出了一種新的脈動陣列結(jié)構(gòu)，但同時為DOA估計的實(shí)時性帶來了障礙，硬件耗時在毫秒級。文獻(xiàn)[11]針對特征分解的繁重計算負(fù)荷以及耗時的問題，對Jacobi算法進(jìn)行了優(yōu)化，提出了一種并行Jacobi算法的一次旋轉(zhuǎn)加速方法，但其硬件耗時仍然在幾百微秒以上。文獻(xiàn)[12]研究了幾種矩陣分解算法在FPGA實(shí)現(xiàn)時對DOA估計性能的影響，其幾種方法的硬件耗時均在微秒級，但沒有考慮完整的MUSIC算法的實(shí)現(xiàn)。到目前對MUSIC算法的FPGA實(shí)現(xiàn)的研究均是基于傳統(tǒng)標(biāo)量陣列，對于基于極化敏感陣列的MUSIC算法的FPGA實(shí)現(xiàn)的研究仍需探索。對DOA估計的實(shí)時性的提高的研究也具有工程應(yīng)用價值。

為了提高硬件效率和算法的實(shí)時性，本文提出一種實(shí)數(shù)化預(yù)處理方法，在此基礎(chǔ)上提出一種基于極化敏感陣列的2維MUSIC算法的FPGA高效硬件實(shí)現(xiàn)方案。從陣列的特點(diǎn)出發(fā)，充分分析陣列的導(dǎo)向矢量，根據(jù)分布式極化敏感陣列的中心對稱結(jié)構(gòu)，得到導(dǎo)向矢量實(shí)部和虛部的對應(yīng)關(guān)系。根據(jù)該對應(yīng)關(guān)系構(gòu)建線性變換矩陣，對接收數(shù)據(jù)做線性變換，使得極化MUSIC算法的計算得到了簡化。FPGA實(shí)現(xiàn)方案主要包含實(shí)數(shù)化預(yù)處理、協(xié)方差矩陣計算、特征值分解、確定噪聲子空間、譜峰搜索幾個模塊。并根據(jù)各模塊獨(dú)立計算的特點(diǎn)設(shè)計了流水線工作模式，充分發(fā)揮了FPGA并行計算的特點(diǎn)，完成了算法的低延遲和資源高效的FPGA實(shí)現(xiàn)。

2 基于極化敏感陣列的MUSIC算法

如圖1所示的陣列為分布式極化敏感陣列，即天線陣列各陣元由1維分量構(gòu)成，且具有不同的極化方式和指向，可無損耗接收對應(yīng)極化方式的來波信號。

圖1 分布式極化敏感陣列示意圖

考慮由N(N=2K或N=2K+1)個陣元組成的一個分布式極化敏感陣列和M個遠(yuǎn)場窄帶完全極化信號，第m個信號的方位角和俯仰角分別為θm和?m， 極化輔助角和極化相位角分別為γm和ηm。假設(shè)入射信號互不相關(guān)且與噪聲之間相互獨(dú)立，則該極化敏感陣列的接收信號可以表示為

式(1)中aθm,?m,γm,ηm為第m個入射信號的導(dǎo)向矢量，任一入射信號的導(dǎo)向矢量具有如式(2)的形式

其中，B=[B1,B2,...,BN]T為極化敏感矩陣，對于不同的極化敏感陣列結(jié)構(gòu)和極化方式，B矩陣具有不同的形式，對于分布式極化敏感陣列，假設(shè)第n個陣元的指向角為αn， 則Bn=[cosαnsinαn]T,n=1,2,...,N。對于第n個陣元空域?qū)蚴噶縰n(θ,?)具有如式(3)的形式

其中，f為信號頻率，τn為不同陣元間的空間延遲，對于平面陣而言，若以原點(diǎn)為參考點(diǎn)，其具體形式為

其中， c=3×108m/s 為光速，(xn,yn)為第n個陣元相對原點(diǎn)的平面坐標(biāo)。

接收信號的協(xié)方差矩陣為

由于信號和噪聲是相互獨(dú)立的，信號的協(xié)方差矩陣可以分解為信號子空間和噪聲子空間兩部分，其中RS是信號的協(xié)方差矩陣，ARSAH是信號部分。對協(xié)方差矩陣進(jìn)行特征分解有

其中，US是 由R的M個較大特征值對應(yīng)的特征向量組成的空間，稱為信號子空間，UN是由R的N－M個較小特征值對應(yīng)的特征向量組成的空間，稱為噪聲子空間[13]。

根據(jù)傳統(tǒng)MUSIC算法可知，信號子空間US與噪聲子空間UN是正交的，即信號子空間的導(dǎo)向矢量也與噪聲子空間正交，那么可以構(gòu)造空間-極化域譜函數(shù)為

式(7)有4個未知變量，那么就需要進(jìn)行一個4維的譜峰搜索才能得到信源的DOA參數(shù)估計，這種搜索有非常多的搜索節(jié)點(diǎn)，在硬件實(shí)現(xiàn)時會造成較高的時間消耗。如果要減少譜峰搜索的搜索時間，那么就需要對極化MUSIC算法譜函數(shù)進(jìn)行降維。

根據(jù)式(2)，式(7)可以進(jìn)一步表示為

其中，hγm,ηm稱 為極化矢量且滿秩，是與DOA參數(shù)無關(guān)的向量，由式(9)表示

因此，極化MUSIC算法的空間譜函數(shù)可以表示為

此過程稱為極化MUSIC算法譜函數(shù)的降維[14]。

3 極化MUSIC算法的FPGA實(shí)現(xiàn)

采用如圖1所示的均勻分布的具有中心對稱結(jié)構(gòu)的分布式極化敏感陣列，在賽靈思(Xilinx)公司的Virtex7-690T FPGA芯片上，使用Verilog或System Verilog語言，在Vivado平臺上進(jìn)行綜合仿真實(shí)現(xiàn)。

3.1 實(shí)數(shù)化預(yù)處理

將省略極化矢量的導(dǎo)向矢量展開為如式(11)的形式

根據(jù)陣元的中心對稱關(guān)系，當(dāng)N為偶數(shù)時，DH的第1個和第K+1個元素分別為

即DH(1)與DH(K+1)的 實(shí)部相反，虛部相同，且DV(1)與DV(K+1)也存在同樣的對稱關(guān)系，對于其他陣元也有同樣的對稱關(guān)系。當(dāng)N為奇數(shù)時，陣列中心對稱的陣元同樣滿足上述對稱關(guān)系。

因此構(gòu)造線性變換矩陣U

其中，IK為K階單位陣，O為K維元素均為0的列向量。對接收信號X(t)做線性變換

此時，Y(t)仍然是復(fù)數(shù)數(shù)據(jù)，但數(shù)據(jù)的實(shí)部已經(jīng)包含了相位信息，另外，硬件實(shí)現(xiàn)時實(shí)部和虛部作為IQ數(shù)據(jù)流分別輸入到系統(tǒng)中，因此，下面本文以Y(t)的實(shí)部為例。圖2顯示了當(dāng)N為偶數(shù)時實(shí)數(shù)化處理的FPGA實(shí)現(xiàn)示意圖。對距離為K的所有元素對，實(shí)部相減，虛部相加，將復(fù)數(shù)數(shù)據(jù)轉(zhuǎn)換為實(shí)數(shù)數(shù)據(jù)。

圖2 實(shí)數(shù)化處理實(shí)現(xiàn)示意圖

而此時導(dǎo)向矢量變?yōu)?/p>

至此，極化MSUIC算法的后續(xù)計算可以完全通過實(shí)數(shù)完成，相比于復(fù)數(shù)計算可以節(jié)約大量的硬件資源，同時縮短了計算時間，為DOA估計的實(shí)時性帶來了可能。

3.2 協(xié)方差矩陣的計算

信號協(xié)方差矩陣R的計算公式為

其中，X為陣列接收信號，L為快拍數(shù)。由于在經(jīng)過實(shí)數(shù)化預(yù)處理后R為實(shí)對稱矩陣，因此只需計算矩陣R的上三角元素的值。此外由于取平均運(yùn)算對后續(xù)計算沒有影響，因此計算時可以省略快拍數(shù)L。綜上可以將協(xié)方差矩陣計算式表示成如式(18)的形式

其中，xl(l=1,2,...,L)為陣列單快拍接收數(shù)據(jù)經(jīng)過實(shí)數(shù)化預(yù)處理得到的數(shù)據(jù)。xlxlT形如式(19)

其中，xln(l=1,2,...,L,n=1,2,...,N)為xl的第n個元素。矩陣共有 ( 1+N)×N/2個上三角元素，每個元素通過定點(diǎn)數(shù)乘法器(Multiplier IP核)及定點(diǎn)數(shù)加法累加器(Accumulator IP核)并行計算，即可得到最終結(jié)果。

3.3 特征值分解

在得到協(xié)方差矩陣后，進(jìn)行相應(yīng)的特征值分解，從而得到特征值和特征向量。Jacobi算法是常用的特征值分解算法之一，主要通過對矩陣做一系列的旋轉(zhuǎn)變換，將其變換成對角矩陣，從而計算出矩陣的特征值和特征向量[15]。本文選用并行Jacobi算法，充分利用FPGA和Jacobi算法的并行計算優(yōu)勢。使用Brent等人[16]提出的脈動陣列結(jié)構(gòu)，將N×N維協(xié)方差矩陣的上三角元素拆分為N-1組，每組為N/2個 2 ×2的矩陣處理單元，根據(jù)每個單元在協(xié)方差矩陣中所處的位置可以分為對角單元和非對角單元，通過并行分組及調(diào)度規(guī)則保證同組N/2個處理單元之間獨(dú)立并行地進(jìn)行Jacobi旋轉(zhuǎn)。當(dāng)陣元數(shù)N=10時，并行Jacobi分組如表1所示。

表1 10×10矩陣的并行Jacobi分組

首先通過對角單元根據(jù)?=1/2arc tan(2b/(a-c))計算旋轉(zhuǎn)角?的值，并根據(jù)式(20)對對角單元進(jìn)行Jacobi旋轉(zhuǎn)運(yùn)算

同時根據(jù)式(21)對非對角單元進(jìn)行Jacobi旋轉(zhuǎn)運(yùn)算

每組Jacobi旋轉(zhuǎn)的FPGA實(shí)現(xiàn)結(jié)構(gòu)如圖3所示。通過調(diào)用CORDIC IP核的arc tan模式計算反正切值，

圖3 單組Jacobi旋轉(zhuǎn)實(shí)現(xiàn)示意圖

從而求得旋轉(zhuǎn)角?，通過調(diào)用CORDIC IP核的sin and cos模式計算旋轉(zhuǎn)角的正弦值和余弦值。完成一組Jacobi旋轉(zhuǎn)之后，需要根據(jù)Jacobi分組進(jìn)行數(shù)據(jù)交換，然后重復(fù)進(jìn)行運(yùn)算，直到完成N－1組Jacobi旋轉(zhuǎn)，即完成1級清掃。由于每一組進(jìn)行相同的迭代計算，為了節(jié)約資源，循環(huán)使用一組硬件資源N－1次來完成清掃，因此需要在輸入端和輸出端進(jìn)行數(shù)據(jù)交換。同樣為了節(jié)約資源，數(shù)據(jù)交換不采用矩陣乘法，而采用狀態(tài)機(jī)直接進(jìn)行寄存器賦值操作。

在進(jìn)行特征值分解時，可根據(jù)實(shí)際需要進(jìn)行多級清掃，清掃的次數(shù)越多特征值矩陣越接近對角陣，即誤差越小，但同時會帶來時間消耗的問題，可根據(jù)實(shí)際需求選擇清掃級數(shù)。本文采用如圖4所示的3級清掃結(jié)構(gòu)，3級清掃的實(shí)現(xiàn)硬件結(jié)構(gòu)相同，為節(jié)省FPGA資源，3級清掃重復(fù)使用同樣的硬件資源。

圖4 特征分解實(shí)現(xiàn)示意圖

3.4 確定噪聲子空間

該模塊首先對特征值進(jìn)行排序，然后根據(jù)排序結(jié)果確定噪聲子空間。對N個特征值進(jìn)行從大到小排序，每個特征值均依次與N個特征值進(jìn)行大小比較，計大于等于此特征值的個數(shù)即此特征值的排序結(jié)果。根據(jù)特征值排序結(jié)果，將對應(yīng)的特征向量構(gòu)成噪聲子空間。由于特征值的物理意義為信號或噪聲的功率，故特征值一定為正數(shù)，故通過運(yùn)算符“<”和“==”進(jìn)行特征值的比較。如圖5為特征值排序的FPGA實(shí)現(xiàn)示意圖。

圖5 特征值排序FPGA實(shí)現(xiàn)示意圖

3.5 2維譜峰搜索

3.5.1 導(dǎo)向矢量計算

實(shí)數(shù)化之后的2維陣列的導(dǎo)向矢量形式為

其中， (θ,?)為搜索角度，根據(jù)搜索范圍及步進(jìn)將其正弦值和余弦值保存到ROM中，導(dǎo)向矢量計算實(shí)現(xiàn)示意圖如圖6所示。

圖6 導(dǎo)向矢量計算示意圖

3.5.2 空間譜函數(shù)計算

由以上分析可知，經(jīng)過降維和實(shí)數(shù)化處理之后的極化MUSIC算法的空間譜函數(shù)可以表示為

其中，U～N為實(shí)數(shù)化預(yù)處理得到的實(shí)數(shù)噪聲子空間。為了簡化計算，構(gòu)造偽譜函數(shù)

計算空間譜函數(shù)值，首先以噪聲子空間的第i列為例計算(θ,?)×N(i)×(i)×Dr(θ,?)，稱為一個內(nèi)部單元，其FPGA實(shí)現(xiàn)結(jié)構(gòu)如圖7所示。經(jīng)過實(shí)數(shù)化處理后，噪聲子空間為N×(N-M)維矩陣，即需要N-M個內(nèi)部單元。如圖8為空間譜函數(shù)計算示意圖。

圖7 內(nèi)部單元計算實(shí)現(xiàn)示意

圖8 空間譜函數(shù)計算實(shí)現(xiàn)示意圖

3.5.3 譜峰搜索

當(dāng)俯仰角和方位角的搜索范圍分別為Rθ和R?，搜索步進(jìn)分別為 ?θ和??時，令I(lǐng)=Rθ/?θ，J=R?/??。 這樣(θ,?) 能夠劃分為I×J的新矩陣，就可以將譜函數(shù)P(θ,?)表 示為I×J的矩陣P

這樣2維極化MUSIC算法的譜峰搜索問題就可以轉(zhuǎn)化為從式(25)表示的2維矩陣的I×J個譜值中搜索局部極小值的問題。

Pi,jP(θi,?j)(i=1,2,...,I,j=1,2,...,J)

若用 表示 ，則需要滿足Pi,j同時小于Pi-1,j-1,Pi-1,j,Pi-1,j+1,Pi,j-1,Pi,j+1,Pi+1,j-1,Pi+1,j,Pi+1,j+1，Pi,j才為一個極小值點(diǎn)，其空間位置譜谷值示意圖如圖9所示。搜索到所有譜谷值位置，然后將其依次按從小到大排序，得到最小的M個譜谷值相對應(yīng)的俯仰角和方位角就是信號DOA參數(shù)的估計值。

圖9 空間譜谷值示意圖

從表達(dá)式可以看出，先將方位角θ固定不變，然后遍歷俯仰角?， 再將θ值進(jìn)行改變，這樣持續(xù)計算，就能夠完成整個空域的搜索。在FPGA實(shí)現(xiàn)時通過FIFO IP核進(jìn)行數(shù)據(jù)傳輸，以實(shí)現(xiàn)圖10所示的九宮格數(shù)據(jù)比較結(jié)構(gòu)。

圖10 數(shù)據(jù)流水線處理框圖

為了進(jìn)一步減少計算量，采用兩級搜索方式。首先進(jìn)行搜索步進(jìn)較大的粗搜索，在整個搜索范圍內(nèi)得到搜索結(jié)果 (,)(m=1,2,...,M)。再將搜索結(jié)果串行輸入到第2級，在粗搜索的結(jié)果適當(dāng)范圍內(nèi)以較小的搜索步進(jìn)進(jìn)行細(xì)搜索，得到最終估計結(jié)果(θm,?m)(m=1,2,...,M)。

3.6 流水線設(shè)計

由于譜峰搜索模塊的耗時大于協(xié)方差矩陣計算、特征值分解和確定噪聲子空間模塊的總耗時，且各個模塊的計算是獨(dú)立的，故采用如圖10所示的流水線處理結(jié)構(gòu)，可以進(jìn)一步提高DOA估計實(shí)時性，但譜峰搜索模塊處理的數(shù)據(jù)為其他模塊上一次工作的輸出，故單次DOA估計時間為譜峰搜索模塊的耗時。

4 實(shí)驗(yàn)結(jié)果與評估

4.1 FPGA資源占用

綜合考慮硬件資源、時間資源以及算法精度之間的平衡，使用定點(diǎn)數(shù)實(shí)現(xiàn)算法，雖然數(shù)據(jù)精度會有一定的犧牲，但在硬件資源、時間資源方面卻可以成倍的節(jié)省。

受FPGA資源的限制，傳統(tǒng)極化MUSIC算法使用8個陣元的分布式極化敏感陣列，而實(shí)數(shù)化的MUSIC算法使用10個陣元的陣列，當(dāng)快拍數(shù)為100，搜索范圍為θ∈[0,360?],?∈[0,30?]，粗搜索和細(xì)搜索的步進(jìn)分別為2°和0.2°，當(dāng)芯片的工作頻率為100 MHz時，消耗的硬件資源情況如表2所示、消耗的時間如表3所示。

表2 資源占用情況

表3 各模塊消耗時間

由表2可以看出，對于查找表(LookUp Table,LUT)、觸發(fā)器(Flip-Flop, FF)和DSP資源，提出的實(shí)數(shù)化的極化MUSIC算法的利用率遠(yuǎn)遠(yuǎn)低于傳統(tǒng)的極化MUSIC算法，這表明提出的方案計算量顯著減少。此外，從表3可以看出，就時間消耗而言，這兩種算法在特征值分解模塊中的差異最大，主要是因?yàn)閭鹘y(tǒng)極化MUSIC算法的復(fù)數(shù)計算使特征值分解中的協(xié)方差矩陣的維數(shù)加倍，而該方案通過實(shí)數(shù)化預(yù)處理避免了這一問題。由于流水線設(shè)計，提出的實(shí)數(shù)化的極化MUSIC算法需要34.60 μs即可輸出一次DOA估計參數(shù)。注意，與傳統(tǒng)的極化MUSIC算法相比，所提出的方案利用了更多的陣元，與基于10陣元陣列的極化MUSIC算法相比，資源占用和時間消耗的改進(jìn)將進(jìn)一步提高。

4.2 算法的DOA估計結(jié)果及性能評估

設(shè)置兩個信源，DOA參數(shù)分別為(77.6?,11.6?)和(277.2?,26.1?)，極化參數(shù)均為(γ,η)=(45?,15?)，如圖11為DOA參數(shù)估計結(jié)果圖。由圖可知，實(shí)數(shù)化的極化MUSIC算法的譜峰要比極化MUSIC算法的譜峰低，但仍然有明顯的譜峰。

圖11 DOA參數(shù)估計結(jié)果

假設(shè)單信源入射到該陣列，信源DOA參數(shù)在搜索范圍內(nèi)隨機(jī)給出，極化參數(shù)為(γ,η)=(45?,15?)，接收信號信噪比為0～28 dB，快拍數(shù)為100，搜索步進(jìn)為0.2°時，使用極化MUSIC算法、通過線性變換實(shí)數(shù)化的極化MUSIC算法以及對接收數(shù)據(jù)線性變換后分別取實(shí)部和虛部的極化MUSIC算法，各做1 000次Monte-Carlo實(shí)驗(yàn)，得到均方根誤差(Root Mean Square Error, RMSE)隨信噪比變化的曲線如圖12所示。其中RMSE的計算式為

圖12 均方根誤差隨信噪比變化曲線

其中，C為Monte-Carlo實(shí)驗(yàn)中測向結(jié)果正確的次數(shù)，θc,?c為第c次實(shí)驗(yàn)的真實(shí)DOA參數(shù)，,為第c次實(shí)驗(yàn)的DOA參數(shù)估計值。由圖12可知，線性變換后的極化MUSIC算法與極化MUSIC算法的性能幾乎一致，取實(shí)部和虛部后性能略有所降低，但隨著信噪比的提高，取實(shí)部后的性能有所提高。相比于實(shí)數(shù)化節(jié)省的計算量，其性能的降低可以忽略，由此得出實(shí)數(shù)化的極化MUSIC算法具有重要的工程應(yīng)用價值。

5 結(jié)論

本文首次研究了極化MUSIC算法在FPGA上的實(shí)現(xiàn)。為了解決傳統(tǒng)復(fù)數(shù)極化MUSIC中的資源占用和時間消耗問題，提出了一種實(shí)數(shù)化預(yù)處理方法，利用均勻圓陣的中心對稱結(jié)構(gòu)，通過線性變換將接收數(shù)據(jù)和導(dǎo)向矢量轉(zhuǎn)換為實(shí)數(shù)。通過實(shí)數(shù)運(yùn)算，協(xié)方差矩陣計算和特征值分解的時間消耗顯著減少。此外，本文還構(gòu)建了所提出的實(shí)數(shù)化的極化MUSIC算法在FPGA上的整體實(shí)現(xiàn)方案，詳細(xì)說明了各模塊的設(shè)計。為了有效地連接所有模塊，本文還設(shè)計了一種流水線實(shí)現(xiàn)結(jié)構(gòu)。根據(jù)實(shí)驗(yàn)結(jié)果，在100 MHz的工作頻率下，完成一次完整的DOA估計需要53.96 μs，同時保證了估計精度。該實(shí)現(xiàn)方案具有精度高、實(shí)時性好、資源消耗低等優(yōu)點(diǎn)，具有較大的工程應(yīng)用價值。