空時自適應處理張量波束成形器的外積合成法

畢權楊，李旦，張建秋

復旦大學 智慧網絡系統研究中心和電子工程系，上海 200433

空時自適應處理(Space-Time Adaptive Processing, STAP)是一種應用于機載雷達的陣列自適應處理技術，可通過雜波協方差矩陣(Clutter Covariance Matrix, CCM)構建自適應的權矢量，來提高信干噪比(Signal to Interference plus Noise Ratio, SINR)，進而提高機載雷達對動目標的檢測性能[1]。在應用中，CCM一般不可知，因此需通過訓練樣本來對其進行估計。研究表明：若訓練樣本滿足獨立同分布(Independent and Identical Distribution, IID)的條件，且樣本數目達到2倍的自由度(Degrees of Freedom, DoF)，則STAP可獲得近似的最優性能[2]。然而由于地表環境以及地形的變化，導致雜波環境往往是非均勻的，因此難以獲得足量的平穩訓練樣本，進而導致STAP的處理性能損失嚴重[2]。

為了解決對足量平穩訓練樣本的要求，文獻報道了許多不同的算法，例如：利用均勻線性陣列(Uniform Linear Array, ULA)下地雜波具有的低秩性，文獻[3-5]分別提出了降秩的STAP(Reduced Rank STAP, RR-STAP)算法，這一算法通過將雜波投影到與雜波正交的子空間，可將訓練樣本的數目降低到2倍CCM的秩。文獻[6-7]分析了RR-STAP算法在不同環境下的性能，指出其在理論上具有很好的應用前景。文獻[8]考慮到傳統RR-STAP算法對數據的矢量化，沒有充分利用多維數據的結構信息，因而將文獻[9]中基于高階奇異值分解(Higher Order Singular Value Decomposition, HOSVD)的子空間估計算法引入到了STAP之中，并提出了一種交替展開高階奇異值分解(Alternative Unfolding HOSVD, AU-HOSVD)算法，從而為STAP算法發展了一種可利用多維數據結構信息的張量處理算法。文獻[10]在文獻[8]的基礎上，完善了AU-HOSVD算法，并將其應用到了多輸入多輸出(Multiple-Input Multiple-Output, MIMO)的STAP算法中。這些基于張量數據處理的算法，發展了傳統的STAP技術，進一步降低了DoF和所需要的訓練樣本數。

本文為文獻[9]中的張量波束成形器推導出一種新的數學表達形式，進而為張量波束成形器的設計提供了一種新的算法：即首先在張量的各個子維度上分別設計出相應的子波束成形器，然后再通過張量外積運算來合成各子波束成形器，以得到STAP所需要的且可充分利用張量數據結構信息的張量波束成形器。分析表明：本文的張量波束成形器與原文獻報道的張量波束成形器相比，在相同DoF基礎上，具有更低的計算復雜度，更有效地去相關處理，使得STAP的性能獲得了進一步的提升，仿真驗證了分析結果的有效性。

本文的組織結構如下：第1節簡單介紹了本文需要的張量基礎以及STAP張量模型；第2節給出了張量波束成形器的外積合成算法；第3節分析和對比了本文提出的張量波束成形器和文獻報道的波束成形器；第4節給出了仿真實驗結果；第5節總結了全文。

1 背景知識

1.1 張量基礎

(1)

式中：?表示等價關系；×n為模n乘積運算符。

式中：vec()為依列將一個多維數組矢量化的運算符;*表示共軛;H表示共軛轉置。

HOSVD是一種常見的張量分解算法，其定義為[11]

(2)

1.2 STAP的張量模型

給定一個均勻線性陣列，假設其天線陣元數為N，且一個相干處理間隔(Coherent Processing Interval, CPI)內的脈沖數為M。稱不同距離門限所接收到的數據為一個快拍，并用x(t)表示第t個距離門限獲得的快拍，其可表示為[13]

(3)

式中：?表示Kronecker積；α(fs,fd)為雜波信號的幅度；ss-t(fs,fd)=ss(fs)?sd(fd)為空時導向矢量；n為均值為零方差為δ2的復高斯觀測噪聲，即n～CN(0,δ2I)，I表示單位陣，并依通常的做法，假設其與雜波在統計上相互獨立；ss(fs)為空間導向矢量；sd(fd)為時間導向矢量，且其分別為

(4)

(5)

(6)

(7)

(8)

2 STAP的張量波束外積合成法

2.1 波束成形器的子波束外積合成法

在陣列信號處理中，波束成形器的目的就是在讓期望信號無失真地通過的同時，盡可能地抑制干擾信號，通常可通過求解下面的優化問題而獲得[17]：

(9)

式中：w為所需波束成形器的權矢量；R=E(x(t)x(t)*)為陣列入射信號x(t)的協方差矩陣；E(·)為期望運算，s0為期望信號的導向矢量。對式(6)的張量模型而言，式(9)的優化問題就變為[18]

(10)

(11)

(12)

將式(12)重寫為張量，則有

(13)

根據式(13)，則可將式(11)重寫為

w1°w2°w3

(14)

2.2 STAP的子波束成形器設計

本節將討論如何設計式(14)中的子波束成形器w1、w2和w3。

(15)

(16)

(17)

(18)

據文獻[20]的降秩理論，對估計得到的CCM做一個低秩近似，就可將DoF降低到雜波子空間的維數。進行這樣的低秩近似后，相應的子波束權矢量就可表示為[13]

(19)

(20)

3 性能分析

本節將對TSS-STAP算法性能進行分析，并將其與最近文獻[10]報道的AU-HOSVD算法進行對比。

3.1 CCM估計

(21)

(22)

將式(21)代入式(22)，有

(23)

(24)

(25)

δ2I

(26)

3.2 算法復雜度

STAP的實時檢測過程，可分為2個部分，第1部分為根據訓練樣本對CCM進行估計，然后通過CCM的特征向量計算投影矩陣，第2部分為通過導向矢量構建自適應的權矢量，然后再利用它們對待檢測快拍進行目標搜索，故本節將主要通過分析計算投影矩陣，以及單次搜索這兩個部分來對比TSS-STAP與AU-HOSVD算法的計算復雜度。

4 仿真實驗

本節將通過仿真實驗來驗證分析結果的有效性，并將與傳統RR-STAP算法以及最新的AU-HOSVD算法進行性能對比。本文將以ULA為例來進行仿真，主要分為兩部分，首先通過數值仿真驗證所提出算法的性能，第1部分為模擬實際環境進行仿真實驗，以驗證算法的有效性。通常，STAP算法的性能用SINR損失因子(SINR Loss Factor)衡量，本文用ρ表示，其定義為輸出SINR與最優STAP算法輸出SINR的比值[10]：

(27)

式中：R為待檢測快拍真實的CCM。

4.1 數值仿真

數值仿真參數設置如下：陣元數N=8，CPI內脈沖數M=24，參數β=1，訓練樣本滿足IID，雜噪比(Clutter to Noise Ratio, CNR)為30 dB，蒙特卡羅仿真1 000次，訓練樣本數P=62，即為RR-STAP理論上所需要的訓練樣本數目(2倍于CCM的秩)。

圖1 CCM的秩隨參數L1的變化曲線

圖3給出了RR-STAP、AU-HOSVD以及模2劃分下TSS-STAP算法的SINR損失曲線圖，其中AU-HOSVD僅給出了最優結果用于對比。如圖3所示，仿真結果中出現了多個主雜波區，其原因為模2劃分下參數β′=L1β，導致了多普勒模糊[23]。此外，圖3也表明：模2劃分下的輸出SINR低于模1劃分，這與前文雜波秩仿真結論一致；本文提出的算法在L1=3時稍優于AU-HOSVD算法，即有稍高的輸出SINR以及更窄的主雜波寬度，當L1=3時，其結果優于L1=2的情況，這與前文的結論模2劃分傾向于更大的L1相一致，但更大的L1將使得多普勒模糊更加嚴重，因此，L1在應用中需要折衷考慮。

圖4為TSS-STAP以及AU-HOSVD算法計算矩陣的運行時間對比，運行環境：CPU為酷睿i7，運行內存為8 G，仿真軟件為MATLAB R2018a。其中圖4(a)為P=62、L2=3時，2種算法計算矩陣的運行時間對比，圖4(b)為P=62、待檢測快拍數為10、參數L2取不同值時，2種算法計算矩陣的運行時間對比。如圖4(a)所示，當L2=3時，本文提出的算法在獲得稍優于AU-HOSVD性能的同時，將運行時間縮短了4～5倍；而根據圖4(b)的結果，可以看到當L2取不同值時，運行時間縮短了4～8倍，這與前文理論分析結果相符。

圖4 計算投影矩陣的運行時間對比

4.2 模擬仿真

下面通過模擬仿真來驗證所提出算法的有效性，主要仿真參數為：陣元數N=8，CPI內的脈沖數M=24，載機飛行高度為1 000 m，飛行速度為225 m/s，脈沖重復頻率為30 kHz，雷達波長為0.03 m，陣元間距為半波長，待檢測地面模擬環境為植被覆蓋環境，天線為正側陣。待檢測目標參數：方位角為40°，歸一化多普勒頻率為0.21，圖5給出了待檢測快拍的功率譜，其中虛線框區域為主雜波區域。

圖5 目標所在快拍的功率譜

在模擬仿真中，AU-HOSVD與TSS-STAP均采用模1劃分，且參數L1=8。圖6給出了P=30時，各算法濾波器的頻率響應，觀察主雜波區域(圖中虛線框部分)對應的頻率響應，可以看出，因為訓練樣本不足，傳統的矢量RR-STAP算法已無法有效抑制主雜波，而AU-HOSVD的雜波抑制沒有TSS-STAP充分。為了更加直觀地對比幾種算法，圖7給出了歸一化多普勒頻率fd以及波達角θ在給定不同值時，各算法的濾波器頻率響應曲線圖，其中圖7(a)對應的fd=0.21，圖7(b)對應的fd=-0.23，圖7(c)對應的θ=41°，圖7(d)對應的θ=-28°，表1則給出了各算法在不同參數下旁瓣增益的平均值。可以發現，當給定的fd與待檢測目標一致時，即fd=0.21，3種算法的旁瓣增益平均值基本相等，濾波器頻率響應基本一致。當fd=-0.23時，RR-STAP的旁瓣高度遠大于另外2種算法，其旁瓣平均增益達到了3.84 dB，而TSS-STAP算法相對于AU-HOSVD算法，雖然主瓣高度相同，但其旁瓣平均增益為-7.65 dB，小于后者的-4.01 dB。而當給定的θ與待檢測目標一致時，即θ=41°，RR-STAP的旁瓣高度稍低于另外兩種算法，其旁瓣平均增益為11.93 dB，當θ=-28°時，根據表1，可知TSS-STAP的旁瓣平均增益最低。

圖6 各算法的濾波器頻率響應

圖7 各算法濾波器頻率響應在給定參數下的對比

表1 各算法旁瓣平均增益對比

圖8為P=30時各算法的檢測結果。可以發現，傳統的矢量RR-STAP算法已無法有效抑制主雜波，而基于張量的低秩濾波器依舊能給出清晰的檢測結果，與前面根據濾波器的頻率響應分析所得出的結果相一致，該結果也驗證了基于張量的STAP方法所需要的訓練樣本數低于傳統的矢量STAP方法，使得其在非均勻雜波應用環境下擁有更好的性能。觀察圖8(b)與圖8(c)，可以發現TSS-STAP相對于AU-HOSVD的濾波結果偽影更少，其原因為：當干擾信號與期望信號的參數不一致時，AU-HOSVD的濾波器頻率響應的旁瓣高于TSS-STAP算法的旁瓣，該結果再次驗證了本文提出的算法能夠更加充分地利用自相關矩陣的全部信息去實現去相關處理。

圖8 各算法的檢測結果

為了定量地描述濾波結果的質量好壞，本文引入圖像熵(Entropy)的概念來進行評價。在圖像處理領域中，圖像熵的大小反映了一副圖像中平均信息量的多少，其定義為[24]

(28)

式中：Pi表示圖像中灰度值為i的像素所占的比例。例如一副純色圖像，其對應的熵為零，那么對于本文中的濾波結果圖，若其中的偽影越少，則灰度值相等的像素點就應該越多，其對應的熵就越低。據式(28)得到圖8(a)對應的熵為6.3，圖8(b)對應的熵為5.4，而圖8(c)對應的熵為4.6。圖9給出了不同訓練快拍數下，TSS-STAP和AU-HOSVD濾波結果對應的熵，該結果表明：本文算法給出的濾波結果中包含了更少的偽影。

圖10為模擬仿真過程中TSS-STAP以及AU-HOSVD算法對待檢測快拍的搜索時間對比，運行環境與前文一致。該結果表明：相對于AU-HOSVD，本文提出的算法在獲得性能提升的基礎上，將搜索時間縮短了約3～4倍，即約為L2倍，這與前文的理論分析相符。

圖9 不同訓練快拍數下濾波結果的熵值對比

圖10 待檢測快拍的搜索時間對比

5 結 論

本文提出了一種新的張量波束成形器設計算法——TSS-STAP算法，通過理論分析，可得出以下幾點結論：

1) 本文方法可在較低DoF的子維度上對張量波束成形器進行設計，因此降低了設計所需要的訓練樣本數。

2) 與文獻報道的張量波束成形器算法相比，本文算法在相同的DoF的基礎上，擁有更低的計算復雜度。

3) 本文提出的算法能更加充分地利用全部自相關矩陣的信息去實現去相關處理，從而使算法性能得到了進一步提升。

仿真結果驗證了上述理論分析的結果，最后本文引入了圖像熵這一概念，用以定量的評價濾波結果的質量。

附錄A

據式(12)的張量信號模型，可將接受到的快拍表示為

(A1)

為了表達簡便，式(12)中的s、L1、L2分別用下標1、2、3代替，同時信號的復幅度被納入在導向矢量中。將式(A1)代入式(8)中，且噪聲與雜波回波信號統計獨立，可計算得：

(A2)

(A3)

(A4)

附錄B

據式(25)，有

(B1)

(B2)

利用Kronecker積的性質[12](A?B)(C?D)=(AC)?(BD)，以及Khatri-Rao積的性質[12](A⊙B)H(C⊙D)=(AHC)(BHD)，將式(B2)中的BTB*展開后，有

(B3)

將式(B3)的結果代入式(B2)中即得到式(26)，推導完畢。