極化敏感陣列到達方向估計方法的FPGA實現

劉魯濤 曹 瑩* 鄭 昱

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

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

1 引言

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

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

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

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

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

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

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

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

其中，B=[B1,B2,...,BN]T為極化敏感矩陣，對于不同的極化敏感陣列結構和極化方式，B矩陣具有不同的形式，對于分布式極化敏感陣列，假設第n個陣元的指向角為αn， 則Bn=[cosαnsinαn]T,n=1,2,...,N。對于第n個陣元空域導向矢量un(θ,?)具有如式(3)的形式

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

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

接收信號的協方差矩陣為

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

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

根據傳統MUSIC算法可知，信號子空間US與噪聲子空間UN是正交的，即信號子空間的導向矢量也與噪聲子空間正交，那么可以構造空間-極化域譜函數為

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

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

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

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

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

3 極化MUSIC算法的FPGA實現

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

3.1 實數化預處理

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

根據陣元的中心對稱關系，當N為偶數時，DH的第1個和第K+1個元素分別為

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

因此構造線性變換矩陣U

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

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

圖2 實數化處理實現示意圖

而此時導向矢量變為

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

3.2 協方差矩陣的計算

信號協方差矩陣R的計算公式為

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

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

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

3.3 特征值分解

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

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

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

同時根據式(21)對非對角單元進行Jacobi旋轉運算

每組Jacobi旋轉的FPGA實現結構如圖3所示。通過調用CORDIC IP核的arc tan模式計算反正切值，

圖3 單組Jacobi旋轉實現示意圖

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

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

圖4 特征分解實現示意圖

3.4 確定噪聲子空間

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

圖5 特征值排序FPGA實現示意圖

3.5 2維譜峰搜索

3.5.1 導向矢量計算

實數化之后的2維陣列的導向矢量形式為

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

圖6 導向矢量計算示意圖

3.5.2 空間譜函數計算

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

其中，U～N為實數化預處理得到的實數噪聲子空間。為了簡化計算，構造偽譜函數

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

圖7 內部單元計算實現示意

圖8 空間譜函數計算實現示意圖

3.5.3 譜峰搜索

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

這樣2維極化MUSIC算法的譜峰搜索問題就可以轉化為從式(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才為一個極小值點，其空間位置譜谷值示意圖如圖9所示。搜索到所有譜谷值位置，然后將其依次按從小到大排序，得到最小的M個譜谷值相對應的俯仰角和方位角就是信號DOA參數的估計值。

圖9 空間譜谷值示意圖

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

圖10 數據流水線處理框圖

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

3.6 流水線設計

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

4 實驗結果與評估

4.1 FPGA資源占用

綜合考慮硬件資源、時間資源以及算法精度之間的平衡，使用定點數實現算法，雖然數據精度會有一定的犧牲，但在硬件資源、時間資源方面卻可以成倍的節省。

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

表2 資源占用情況

表3 各模塊消耗時間

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

4.2 算法的DOA估計結果及性能評估

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

圖11 DOA參數估計結果

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

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

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

5 結論

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