基于張量結構的快速三維稀疏貝葉斯學習STAP方法

崔 寧 行 坤 段克清 喻忠軍*

①(中國科學院空天信息創新研究院 北京 100094)

②(中國科學院大學電子電氣與通信工程學院 北京 100049)

③(中山大學電子與通信工程學院 廣州 510006)

1 引言

空時自適應處理(Space-Time Adaptive Processing,STAP)作為一種實用的雜波抑制技術，可有效地改善機載雷達在強雜波環境下的運動目標檢測能力。該技術主要利用雜波方位和多普勒之間的耦合特性，在空時二維平面內構建最優自適應濾波器，因此具有比傳統多普勒濾波器更好的雜波抑制能力。然而，待檢測區域實際雜波分布特性通常是未知的，需要利用待檢測距離單元(Cell Under Test,CUT)相鄰數據作為訓練樣本經似然估計獲取。為分析雜波協方差矩陣(Clutter Covariance Matrix,CCM)估計誤差對STAP性能影響，Reed,Mallet和Brenna共同提出RMB準則，即輸出信雜比相較于最優性能損失不超過3 dB所需最小獨立同分布(Independent and Identically Distributed,IID)訓練樣本數為2倍系統自由度(Degrees Of Freedom,DOF)[1]。

當機載雷達處于非正側工作時，雜波多普勒頻率與空間方位關系隨距離改變而變化，即存在雜波非平穩現象。同時，由于機載預警雷達常工作于中/高脈沖重復頻率(Pulse Repetition Frequency,PRF)模式，因此雷達回波往往存在嚴重的距離模糊。上述非平穩雜波和距離模糊交織在一起，進而導致后續傳統二維STAP (Two-dimensional STAP,2D-STAP)處理存在以下問題：(1)各距離單元雜波分布特性不同，使得用于估計CCM的IID訓練樣本數嚴重不足，引起STAP性能惡化；(2)CUT中存在多個具有不同多普勒頻率的方位主瓣雜波，大大增加了盲速范圍。為消除非平穩雜波影響，國外學者陸續提出各類雜波譜補償方法，如多普勒補償法、高階多普勒補償法、角度-多普勒補償法、自適應角度-多普勒補償法等[2]，通過各類數學變換，將不同距離單元雜波譜校正一致。然而，上述補償類方法無法同時補償由距離模糊導致的重疊雜波，因此僅適用于低PRF工作模式[3]。

三維STAP(Three-dimensional STAP,3DSTAP)在方位-多普勒二維處理的基礎上利用平面陣列增加了俯仰維度自適應處理能力，可在距離模糊情況下取得很好的雜波抑制效果。然而，更多的系統DOF令3D-STAP處理所需IID樣本數遠大于傳統2D-STAP，為解決上述問題，降維技術(Reduced-Dimensional,RD)由二維拓展至三維，形成3D-RDSTAP方法[4]。RD處理可將IID樣本需求降至2倍降維后的系統DOF，但在極端非均勻雜波環境下依然難以獲得充足的訓練樣本。近年來，與稀疏恢復(Sparse Recovery,SR)技術相結合的STAP處理方法受到了廣泛關注，大量研究表明SR-STAP只需極少訓練樣本甚至單樣本即可獲得良好的雜波抑制能力[5,6]。最近有學者將SR應用于3D-STAP并提出3D-SR-STAP方法[7]，結果表明僅需數個樣本3DSR-STAP即可實現非平穩雜波的有效抑制。雖然SR-STAP在小樣本條件下相比傳統方法具有顯著的優勢，但巨大的計算開銷嚴重制約了它的實際應用。目前，已有學者針對SR-STAP計算開銷問題提出許多快速算法，例如，結合降維方法形成局域稀疏處理[8]，通過譜先驗知識輔助對空時格點預篩選[9]以及設計SR-STAP快速處理算法[10]。然而，以上方法均是對矢量化后的回波數據進行計算，所需處理的字典維度將隨各維度數據增加而急劇上升，對硬件處理平臺的計算負載和存儲能力提出了極高要求[11]。

張量作為一種處理多維數據的高效方式，已在機器學習[12]和圖像處理[13]等領域取得了廣泛的研究和應用。不同于傳統矢量基處理，張量處理借助數據Kronecker結構特性，將大規模矩陣運算拆解為多個小尺度矩陣運算，避免因重復計算而產生的額外開銷。據此，本文結合平面陣機載雷達回波三階結構特性提出一種基于張量的快速三維稀疏貝葉斯學習STAP (3D Sparse Bayesian Learning STAP,3D-SBL-STAP)方法。首先，建立基本的三維信號模型；然后在此模型基礎上，本文將2D-SR-STAP問題拓展至3D-SR-STAP問題，并給出了多觀測樣本矢量(Multiple Measurement Vectors,MMV)下3D-SR-STAP優化模型；接著對所提張量MSBL方法的數學原理進行描述，并從運算復雜度對所提方法的性能優勢進行分析；最后，數值實驗驗證了所提方法。

2 三維信號模型

如圖1所示，載機平臺在高度H處以恒定速度vp沿y軸正向飛行，其上裝配M ×N個理想全向陣元構成的均勻平面機載相控陣雷達，任意相鄰陣元間距d均為半波長。假定待探測區域位于天線遠場，則對于不同陣元，其來向回波可視為具有相同入射角。在每個相干處理間隔內，雷達通過收/發天線以恒定的PRF向外發射K個窄帶調頻相參脈沖串并接收回波。θ和φ分別表示陣列與地面雜波塊之間的方位角和俯仰角。θp為天線陣列擺放方向與載機飛行方向之間的夾角。

圖1 平面陣機載雷達幾何照射圖Fig.1 The geometry of planar phased-array airborne radar

經T/R組件和A/D采集處理后，NMK×1基帶回波信號矢量x可表示為

其中，xc代表雜波成分；Nr和Nc分別表示距離模糊數和單個距離環內所劃分的雜波塊數；ξi,j為第i個模糊距離的第j個雜波塊所對應的復反射系數；n通常假定為零均值、方差σ2的復高斯白噪聲；s(θ,φ)∈CNMK×1為三維空時導向矢量，可進一步展開寫為

其中

其中，se(φ),sa(θ,φ)和sd(θ,φ)分別為俯仰、方位和多普勒導向矢量；fr代表雷達PRF；λ為雷達工作波長； (·)T為轉置操作；?表示Kronecker積。進而CCM可表示為

其中，E{·},(·)?,|·|和 (·)H分別為期望、共軛、取模和共軛轉置操作。實際中，由于先驗雜波分布未知，常使用CUT相鄰數據作為訓練樣本進行最大似然估計(Maximum Likelihood Estimation,MLE)來替代Rc，即,L為訓練樣本數。最優三維空時自適應處理優化問題定義為

其中，st代表目標方向的三維空時導向矢量。求解上述問題得最優3D-STAP權值為

其中，(·)?1為矩陣求逆操作。當MLE樣本不足時，無 法直接求逆，此時可通過對角加載改善性能[14]。

3 三維稀疏STAP基本原理

首先，對方位-俯仰-多普勒三維連續空間做網格離散化處理，將其均勻劃分為Na×Me×Kd個空間格點。其中，NaβaN,MeβeM和KdβdK分別代表方位、俯仰和多普勒維劃分網格數；βa,βe和βd為各維度的網格劃分系數。所有格點對應的三維空時導向矢量可組成過完備字典S，具體可寫為

其中，ξ[ξ1ξ2...ξNaMeKd]T為反射系數矢量。由于在方位-俯仰-多普勒三維空間內雜波呈現稀疏分布特性，可知ξ中大部分元素為0[5]，滿足稀疏假設，因此可通過SR來估計雜波區域。3D-SRSTAP可表述為如下優化問題

其中，diag(·)表示將矢量排成對角矩陣；δ為對角加載值；I為單位對角矩陣。在獲得估計CCM后，可進一步通過式(5)計算STAP權值。上述模型基于單觀測矢量(Single Measurement Vector,SMV)，但實際中常采用MMV模型以進一步減小噪聲影響[15]。MMV模型下的3D-SR-STAP優化問題可描述為

其中，‖·‖2,1是l2/l1混合范數，具體為對空時維度(列)使用l1范數，對距離維度(行)使用l2范數；為斐波那契范數，用來度量恢復矩陣和原始矩陣之間的相近程度；Ξ[ξ1,ξ2,...,ξL]代表恢復矩陣。相應地，MMV模型下的構造式為

4 張量基處理結構分析

常規矢量基處理先將高維回波X ∈CN×M×K排為NMK×1的向量x，然后再進行處理。矢量處理實現簡單，但在計算過程中，單維字典Sa,Se和Sd通過Kronecker乘積耦合，令完備字典S尺度擴充至NMK×NaMeKd，由此產生的大量冗余項給后續操作帶來嚴重的計算負擔。張量基運算直接在多維數據域處理，無須構建高維字典S，從而避免由此導致的重復計算。具體地，對3階平面陣數據而言，在進行乘法操作時，X沿方位、俯仰和多普勒維分別與單維字典Sa,Se和Sd進行3次小規模矩陣乘法，每次待處理矩陣的維度僅與上一次處理有關，中間無重復項生成，因此計算效率更高。下面將對張量處理過程展開詳細分析。

其中，ZFNa,Me,Kd,L{ξ}為待乘張量；(·)×n為張量的n階模乘。以方位-俯仰-多普勒模乘順序為例，張量模乘的具體實現過程如圖2所示。首先在方位模乘中，將Z按方位模作切片展開，即ZMN,MK{Z}，然后作矩陣相乘得Z′，所需乘法次數為NaNMK；依次在俯仰模乘和多普勒模乘中采取相同操作，所需乘法次數分別為MeMNaK和KdKNaMe；最后可得3階張量模乘總的乘法次數為NaK(NM+MeM+KdMe)。

圖2 張量基處理流程Fig.2 The flow of tensor-based processing

通過上述分析，可知張量模乘將矢量基乘法SHξ拆解為3個小尺度矩陣相乘：和從而令所需乘法次數NMKNaMeKd中部分乘積項變為求和項NaK(NM+MeM+KdMe)，進而使計算復雜度由O(NMKNaMeKd)降至O(NaKKdMe)。進一步如圖3所示，對比兩種處理方式可發現在矢量基乘法中存在大量由Kronecker乘積產生的冗余項，而張量運算由于分維度實現，不產生額外的冗余項，也就避免無用的運算開銷。

圖3 兩種不同結構對比Fig.3 The comparison of two different structures

5 三維張量稀疏貝葉斯學習STAP

假定Ξ中元素相互獨立，且每列均服從零均值復高斯先驗分布

結合式(13)和式(14)，由貝葉斯準則可得后驗概率pΞ|X,Γ,σ2也服從復高斯分布CN(Y,Σ)[16]，其中

進而Ξ可由均值Y估計。未知參數采用期望最大化方法計算，求解過程如下：首先，基于對數最大似然構造損耗函數

其中，Cσ2I+SΓ SH。利用矩陣恒等變換，式(16)可改寫為

然后，分別令其對γ和σ2的1階偏導數為0，得每次迭代的最優估計為

其中，DiΣi,i/γi;Tr(·)代表矩陣的跡。以上為MSBL方法，下面引入張量形式。γi估計過程僅使用Σi,i，因此將其放入循環內部求解

其中，QΓ SHC?1,TSΓ。經簡化處理后，復雜度最高步驟僅余矩陣C和Y的計算。首先，考慮C的求解，將T排為4階張量T則

其中，locd,loca和loce分別代表第i個三維空時導向矢量對應字典Sd,Sa和Se中的列。式(20)等效矩陣實現方式如下：

由于采用矩陣存儲，式(21)具有更快的求解速度。此外，建議按照維度由小到大順序計算以獲得最高效率。類似的，再次整理4階張量T為NMK×Na×Me×Kd，則T SH的張量形式為

至此矩陣C張量重構完畢，下面將闡述Y的求解。

由于Γ SH不具備Kronecker結構無法直接分解，因此調整矩陣運算順序為先計算SHC?1，再計算Γ與其相乘。將C?1表示為C?1∈CN×M×K×NMK，進而Y計算可改寫為

其中，⊙為Hamdamard乘積；為系數張量，定義如下：將γ復制NMK次生成NaMeKd×NMK系數矩陣D，緊接著把D折疊為4階系數張量D。張量MSBL (Tensor MSBL,TMSBL)整體流程如表1所示，其中(·)i?1代表前一次迭代結果。進一步歸納可知張量處理前提為待相乘的兩個矩陣至少一個具備Kronecker結構，滿足使用條件即可采用張量處理。

表1 TMSBL算法Tab.1 TMSBL algorithm

6 計算復雜度分析

本節將TMSBL算法與多種SR方法的計算復雜度進行比較，其中正交匹配追蹤(Orthogonal Matching Pursuit,OMP)[17]為貪婪類算法，具有極快計算效率；多迭代自適應(Multiple Iterative Adaptive Approach,MIAA)[18]和多欠定系統聚焦式求解(Multiple Focal Underdetermined System Solver，MFOCUSS)[19]為SR常用方法；MSBL[16]為較新方法，其不依賴參數，因此具有較高的穩健性；多樣本快速收斂SBL (Multiple Fast Converging SBL,MFCSBL)[10]是在MSBL基礎上改進，因此具有更快求解速度。以L個訓練樣本求解所需復乘法次數作為計算復雜度量化統計，同時由于后續構造CCM步驟相同，這里為簡化分析不計入此步操作。表2為各算法計算復雜度，其中OMP算法稀疏度以平面陣雜波秩rsmin(M,Nse)rsa來近似估計[4]，rsa和Nse分別為方位線陣雜波秩和俯仰維子陣數。KOMP,KMIAA,KMFOC,KMSBL,KMFCSBL和KTMSBL分別代表各算法單樣本求解所需迭代次數，通常MFOCUSS和MSBL迭代次數要大于MIAA。此外，由于張量處理不影響算法固有性能，可知KMSBLKTMSBL。一般而言SR-STAP僅需數個樣本即可有效恢復，因此L要遠小于其他參數。

表2 計算復雜度Tab.2 Computational complexity

從表2的統計結果來看，OMP由于只涉及低復雜度的矩陣乘法，具有最高的計算效率；MIAA,MFOCUSS 和MSBL 計算復雜度中分別包含NaMeKd的一次項、二次項和三次項，因此運算開銷呈指數級依次增加。雖然OMP求解的復雜度要遠低于后幾種SR算法，但局部貪婪式搜索性能在非正側視情況下會急劇惡化，無法滿足應用要求；而后幾種方法盡管估計精度較高，但高度復雜的求解過程對實際硬件處理提出了很高的要求，難以實現。所提TMSBL在運算效率上復雜度量級為略高于OMP的O(NMKNaMeKdLKOMP)，這表明TMSBL所需求解時間接近于OMP的量級；在性能上，TMSBL通過張量結構改進MSBL方法，之前分析已表明張量處理加速核心是利用Kronecker性質來避免冗余項計算，所以TMSBL具有與MSBL相同的處理性能。綜上分析可知TMSBL借助張量計算完美地解決了速度和性能之間的矛盾，因此相比其他SR算法具有更顯著優勢。

7 數值實驗

實驗基于仿真回波數據來驗證所提3 DTMSBL-STAP方法的性能。在仿真過程中，天線陣列工作在前視模式，具體仿真參數如表3所示。代碼中涉及的張量計算由MATLAB內部Tensor工具箱實現[20]。

表3 雷達系統參數Tab.3 Radar system parameters

7.1 三維譜估計結果比較

本節檢驗TMSBL的雜波譜估計效果，其中三維離散網格劃分選取βa4,βe2和βd4，迭代終止門限μ10?5，訓練樣本數L6。以CCM的最小方差無失真響應(Minimum Variance Distortionless Response,MVDR)譜作為對比，定義如下：

其中，s(fa,fe,fd)為對應頻點的三維空時導向矢量，P(fa,fe,fd)為譜強度。

圖4為確知CCM的MVDR譜和TMSBL稀疏恢復譜。由圖4(a)可知待估計雜波主要包含兩部分：近程雜波和交疊的遠程雜波。從圖4(b)TMSBL的估計結果來看，在小樣本條件下不同雜波區域均能有效恢復。當繼續增加格點劃分密度至βa8,βe4和βd8時，如圖4(c)所示，由網格失配導致的雜波 展寬區域明顯縮小，從而估計精度得到提升。

圖4 三維雜波譜估計結果Fig.4 The estimation of three-dimensional clutter spectrum

7.2 稀疏方法雜波抑制性能對比

下面將3D-TMSBL-STAP的雜波抑制性能分別與現有3D-OMP-STAP,3D-MIAA-STAP,3DMFOCUSS-STAP,3D-MFCSBL-STAP和3D-MSBL-STAP方法對比，其中網格劃分、迭代條件和樣本選取均與之前相同。特別地，對于3D-OMPSTAP，以雜波秩作為稀疏度估計，設置其略大于rs。雜波抑制性能使用信雜噪比(Signal-to-Clutterand-Noise Ratio,SCNR)損失評價，定義如下[21]：

圖5為信雜噪比損失結果。在非平穩雜波條件下，3D-OMP-STAP處理性能最差，這是由于較大的網格失配導致貪婪類算法在搜索過程中受到嚴重干擾，從而性能損失嚴重。與之對比，3D-SR-STAP方法即使在該條件下依然能取得較好的抑制效果，其中SBL類方法性能最優，信雜噪比損失結果要好于其他SR方法。需要注意的是，所提3D-TMSBLSTAP方法與3D-MSBL-STAP方法性能一致，說明張量處理對算法固有雜波抑制性能無損失。

圖5 不同SR-STAP方法的SCNR損失結果Fig.5 The SCNR loss results of different SR-STAP methodss

7.3 傳統方法雜波抑制性能對比

本節對比傳統3D-STAP方法與3D-TMSBLSTAP方法在小樣本L6條件下的雜波抑制性能，所選取的3D-STAP方法包括：樣本矩陣求逆(Sample Matrix Inversion,SMI)[1]、局域聯合處理(Joint Domain Localized,JDL)[22]、多級維納濾波(Multistage Wiener Filter,MWF)[23]、直接數據域(Direct Data Domain,DDD)[24]和子陣合成(Subarray Synthesis,SS)[4]。具體地，JDL在方位-俯仰-多普勒向選取5×6×5波束；MWF的分解次數設為rs；DDD的滑動窗口選取各維度的一半；SS合成行列陣元數為M和N的十字子陣。

從圖6可以看到，傳統3D-STAP方法在樣本短缺情況下形成凹口遠寬于實際主雜波區域，此時會對主雜波區域外目標強度造成嚴重損失。雖然3DDDD-STAP方法僅靠單樣本即可構建雜波濾波器，但由于存在孔徑損失，可以看到非主雜波區域性能下降嚴重。而3D-TMSBL-STAP方法基于SR理論，對樣本數依賴并不嚴重，在小樣本條件下具有遠好于其他傳統方法的雜波抑制性能。

圖6 傳統STAP方法的SCNR損失結果Fig.6 The SCNR loss results of conventional STAP methods

如圖7所示，以信雜噪比損失的平均值(平均SCNR損失)作為評價[4]，當可用訓練樣本數增加時，傳統方法性能逐漸改善。JDL和MWF兩種方法由于降維損失在L大于50后開始趨于穩定，此時平均信雜噪比損失相比于3D-TMSBL-STAP方法依然相差4～5 dB。而SS和SMI方法盡管最終性能要好于3D-TMSBL-STAP方法，但需要充足樣本才能實現，不適應于嚴重非均勻雜波環境。以上結果進一步說明3D-TMSBL-STAP方法相比于傳統3D-STAP方法更適用于非均勻雜波環境。

圖7 平均SCNR損失結果Fig.7 The result of average SCNR loss

7.4 運行時間對比

最后將對比不同N,M和K大小下3D-TMSBLSTAP方法與現有3D-SR-STAP算法的實際運行時間。參數設置與之前相同，不再贅述。共選取7組數據，對應維度為俯仰×方位×多普勒，相鄰數據僅在單維尺度上變化。各算法均執行于MATLAB R2019a，硬件配置為Intel Core i9-9900KF 3.6 GHz加64 GB內存，采用MATLAB內部計時器統計運行時間。

如圖8所示，從整體運行時間上看，MSBL要遠高于其他方法1～3個數量級，原因是計算過程中涉及兩個NaMeKd×NaMeKd的大規模矩陣乘法計算。緊接著是MIAA,MFCSBL和MFOCUSS，三者的運行時間處于中等水平。運算開銷最小的是OMP和TMSBL，其中TMSBL運算時間要小于OMP，這是因為當工作于大尺度數據時，OMP受rs約束迭代次數較多。綜上所述，可知TMSBL借助 張量處理取得了顯著的速度提升。

圖8 不同數據大小下實際運行時間Fig.8 The running time of different input data size

8 結論

3D-SR-STAP僅需數個樣本即可實現有效的非平穩雜波抑制，但巨大的計算開銷限制了它的應用。本文從空時字典的Kronecker結構特性出發，提出基于張量處理的3D-MSBL-STAP方法，通過低復雜度張量操作加速求解過程，從而獲得遠快于矢量基方法的運算速度。仿真結果表明在本文條件下，所提方法可在維持原算法性能的同時，將計算速度提升2～3個數量級，因此更具實際工程價值。此外，還可將所提張量處理方法拓展應用于其他矢量基求解算法，以提升其運算速度。未來，本團隊將改進張量的讀寫方式，以提升其處理速度，并嘗試融入并行處理，以期獲得進一步加速。