章濤，鐘倫瓏，來燃，郭駿騁

中國民航大學 天津市智能信號與圖像處理重點實驗室，天津 300300

空時自適應(yīng)處理(Space-Time Adaptive Processing, STAP)是機載陣列雷達抑制雜波的一種有效方法[1-2]。STAP的雜波抑制性能主要取決于雜波協(xié)方差矩陣(Clutter-plus-noise Covariance Matrix, CCM)的估計精度。傳統(tǒng)STAP方法采用統(tǒng)計估計的方法獲得CCM估計值，若需保證輸出信雜比相較最優(yōu)值下降不超過3 dB，則至少需要系統(tǒng)自由度2倍數(shù)量的獨立同分布(Independent Identically Distributed, IID)雜波樣本[3]。然而在實際系統(tǒng)中，當雜波呈現(xiàn)非平穩(wěn)或非均勻特性時，很難獲得足夠的IID樣本，進而使得STAP雜波抑制性能嚴重損失[4]。

為了減少STAP所需雜波樣本數(shù)量，研究人員提出了多種方法，主要包括降維STAP方法和降秩STAP方法。降維STAP方法利用與數(shù)據(jù)無關(guān)的變換降低接收數(shù)據(jù)維數(shù)，使得所需雜波樣本數(shù)量為降維后數(shù)據(jù)維數(shù)的2倍，主要方法包括：輔助通道處理(Auxiliary Channel Processing, ACP)方法[5]、擴展因子化(Extended Factored Approach, EFA)方法[6]、廣義多波束(Generalized Multiple Beams, GMB)方法[7]等，但這類方法不可避免地存在系統(tǒng)自由度損失。降秩STAP方法則利用與數(shù)據(jù)有關(guān)的變換對接收數(shù)據(jù)進行處理，使得所需雜波樣本數(shù)降低為雜波秩的2倍，主要方法包括：正交投影處理(Orthogonal Projection Processor, OPP)方法[3]、最小功率特征對消(Minimum Power Eigen canceller, MPE)方法[8]、多級維納濾波(Multistage Winer Filter, MWF)方法[9]等，但是這類方法的性能依賴于雜波秩，在實際雜波環(huán)境中很難準確估計雜波秩。

近年來，稀疏恢復(fù)方法在信號處理領(lǐng)域快速發(fā)展[10-11]。由于稀疏恢復(fù)方法可以使用很少樣本精確地恢復(fù)未知稀疏信號，已成為機載雷達雜波抑制的研究熱點[12-13]。目前比較有代表性的雜波空時譜稀疏恢復(fù)方法包括：文獻[14]提出了利用1范數(shù)最小化實現(xiàn)空時譜稀疏恢復(fù)的凸優(yōu)化方法。文獻[15-17]提出了基于p范數(shù)(0

稀疏恢復(fù)STAP方法將雜波信號看作在某些已知基向量下是稀疏的，這些基向量被稱為稀疏恢復(fù)字典。由于空時平面被離散地劃分為若干個網(wǎng)格點來構(gòu)建空時導(dǎo)向矢量字典，離散化的字典與連續(xù)的雜波譜參數(shù)間的誤差問題被稱作字典失配問題，嚴重影響稀疏恢復(fù)效果。在正側(cè)視陣情況下，若稀疏恢復(fù)空時平面網(wǎng)格劃分與雜波脊斜率不匹配；或是當陣列為非正側(cè)視架設(shè)，雜波脊呈曲線形式時，則雜波大多不會位于預(yù)先離散化的空時平面網(wǎng)格點上，字典失配將嚴重影響稀疏恢復(fù)STAP的雜波抑制性能。雖然增加空時導(dǎo)向矢量劃分密度可以減少字典失配的影響[19]，但字典網(wǎng)格間隔越小，空時導(dǎo)向矢量的相關(guān)性越強，會導(dǎo)致稀疏恢復(fù)性能下降，而且計算量也將大幅增加[10]。

目前，針對稀疏恢復(fù)中字典失配問題的研究主要集中在一維離散參數(shù)字典的補償方面，如：文獻[20]提出了一種稀疏總體最小二乘(Sparse Total Least-Squares, STLS)方法，結(jié)合動態(tài)字典模型解決一維波達方向(Direction of Arrival, DOA)估計中的字典失配問題。文獻[21]提出利用基追蹤方法解決DOA估計中的字典失配問題。為了避免正則化參數(shù)難以確定的問題，文獻[22] 提出了一種利用動態(tài)字典解決字典失配問題的稀疏貝葉斯學習(Off-Gird Sparse Bayesian Inference, OGSBI)方法，但該方法僅能處理一維離散參數(shù)字典失配問題，不能直接用于存在二維字典失配的雜波空時譜稀疏恢復(fù)。目前針對雜波空時譜稀疏恢復(fù)中二維離散參數(shù)字典失配問題的研究還比較少，且有一定的適用局限性。例如，文獻[23]在2013年首次提出了雜波空時譜估計中存在字典失配問題，并給出了一種利用字典學習的失配補償方法，但該方法僅適用于正側(cè)視陣情況。文獻[24]提出了利用載機平臺運動參數(shù)獲得雜波脊先驗信息后細化字典間隔來解決字典失配問題的方法，但這種方法依賴載機平臺傳感器，有些情況并不適用。文獻[25]提出了基于正交匹配追蹤(Orthogonal Matching Pursuit, OMP)的方法，通過梯度下降方法尋找與真實模型匹配的字典向量來解決字典失配問題，但OMP方法的性能對參數(shù)選擇依賴性較高。

針對目前已有的字典失配情況下稀疏恢復(fù)STAP算法存在僅適用于正側(cè)視陣、依賴機載平臺其他傳感器、算法性能對參數(shù)選擇敏感等問題，本文將OGSBI方法擴展到二維字典失配的情況，提出了一種解決空時二維字典失配的稀疏貝葉斯學習STAP方法。本方法的主要創(chuàng)新性包括：① 針 對空時二維參數(shù)字典失配，建立了一種利用二維泰勒級數(shù)的空時二維字典失配誤差動態(tài)模型；② 將空時二維字典失配誤差作為超參數(shù)，給出了期望最大化(Expectation-Maximization, EM)方法迭代估計公式，進而修正空時二維參數(shù)字典失配誤差, 最后利用修正后的空時導(dǎo)向矢量字典估計雜波協(xié)方差矩陣及雜波空時譜。實驗表明，本文方法能夠有效提高雜波譜稀疏恢復(fù)精度，STAP處理性能優(yōu)于已有字典預(yù)先離散化處理的稀疏貝葉斯學習STAP算法。

1 信號模型

考慮采用均勻線陣的機載陣列雷達，如圖1所示，陣列天線由M個陣元組成，陣元間距d=λ/2，λ為雷達工作波長；H為載機平臺高度；vp為速度，且沿x軸運動；α為雜波散射點P與飛行方向間的夾角；θ和γ分別為俯仰角和方位角；φ為陣列軸線與飛行方向間的夾角，即當φ=0°時為正側(cè)視陣模式，φ=90°為正前視陣模式。

圖1 機載雷達陣列幾何結(jié)構(gòu)

雷達在每個相干處理周期內(nèi)發(fā)射N個脈沖，脈沖重復(fù)頻率為fr。當不存在距離模糊時，對于不存在目標的距離單元的回波進行采樣可以獲得一個雜波觀測矢量數(shù)據(jù)(通常稱為一個空時快拍數(shù)據(jù))x∈CNM×1

x=xc+n

(1)

式中：xc為雜波分量；n為噪聲分量。每個距離單元的雜波可以看成是多個雜波散射體空時響應(yīng)的疊加，即

(2)

式中：Nc表示雜波散射體數(shù)量；ai、fd,i和fs,i分別表示復(fù)幅度、多普勒頻率和空間頻率。s(fd,i,fs,i)∈CNM×1稱為空時導(dǎo)向矢量，即s(fd,i,fs,i)=sd,i(fd,i)?ss,i(fs,i)，?表示Kronecker積運算。其中:

空間導(dǎo)向矢量ss,i(fs,i)定義為

ss,i(fs,i)=[1，exp(j2πfs,i)，…，

exp(j2π(M-1)fs,i)]T

(3)

時間導(dǎo)向矢量sd,i(fd,i)定義為

sd,i(fd,i)=[1,exp(j2πfd,i),…,

exp(j2π(N-1)fd,i)]T

(4)

空間頻率fs,i滿足

fs,i=cos(γi-φ)

(5)

多普勒頻率fd,i滿足

(6)

在歸一化空間頻率和多普勒頻率構(gòu)成的空時平面上，雜波脊的斜率(折疊系數(shù))為

(7)

CCM定義為

R=E(xxH)

(8)

在噪聲向量n假設(shè)為零均值高斯分布的情況下，基于最大信雜噪比(Signal-to-Clutter-plus-Noise Ratio, SCNR)準則，可以得到自適應(yīng)權(quán)矢量的STAP輸出為

y=wHx

(9)

式中：自適應(yīng)權(quán)矢量w可以通過CCM計算獲得，即

(10)

其中：st(fd,fs)為待檢測單元的空時導(dǎo)向矢量。由于CCM未知，實際系統(tǒng)中可以通過與目標檢測單元雜波分布相同的雜波樣本計算CCM的估計值，如

(11)

其中：訓練樣本數(shù)L需大于系統(tǒng)自由度的2倍，這種方法稱作采樣矩陣求逆(Sample Matrix Inversion, SMI) STAP方法[1]。但對于非均勻、非平穩(wěn)的雜波環(huán)境不同距離單元雜波數(shù)據(jù)難以構(gòu)成足夠的IID樣本，因此這種方法在有些情況下并不適用[12-13]。

2 空時功率譜稀疏恢復(fù)及字典失配問題

2.1 空時功率譜稀疏恢復(fù)

基于雜波譜稀疏恢復(fù)的STAP方法將上述方法中的CCM估計通過稀疏恢復(fù)的方法實現(xiàn)，可有效減少所需IID雜波訓練樣本數(shù)[12-13]。這類方法首先將雜波譜空時平面劃分為Ns×Nd的網(wǎng)格。其中，Ns=ρsM,Nd=ρdN,ρs和ρd為網(wǎng)格劃分系數(shù)，且ρs>1,ρd>1。網(wǎng)格點所對應(yīng)的空時導(dǎo)向矢量集合可以表示為

Ψ=[ψ1,ψ2,…,ψNs×Nd]=Sd?Ss

(12)

式(12)在稀疏恢復(fù)STAP中被稱為空時導(dǎo)向矢量字典，其中

Ss=[ss,1,ss,2,…,ss,Ns]∈CM×Ns

(13)

Sd=[sd,1,sd,2,…,sd,Nd]∈CN×Nd

(14)

單快拍(Single Measurement Vector, SMV)雜波樣本數(shù)據(jù)x的稀疏恢復(fù)模型可以表示為

x=Ψa+n

(15)

式中：a=[a1,a2,…,aNsNd]T為稀疏恢復(fù)支撐集向量，其每一個非零元素表示一個對應(yīng)的雜波單元。通過稀疏恢復(fù)理論，支撐集向量a可以通過下面最優(yōu)化方法獲得

(16)

通過稀疏恢復(fù)方法獲得支撐集向量a后，CCM矩陣可通過式(17)重構(gòu)獲得

(17)

為了提高稀疏恢復(fù)性能，可以利用多快拍(Multiple Measurement Vector, MMV)雜波樣本數(shù)據(jù)進行聯(lián)合稀疏恢復(fù)。K個雜波樣本數(shù)據(jù)X=[x1,x2,…,xK]∈CNM×K可以表示為

X=ΨA+N

(18)

式中：A=[a1,a2,…,aK]為稀疏恢復(fù)的支撐集矩陣，其每個非零行表示一個對應(yīng)的雜波單元；N=[n1,n2,…,nK]為噪聲矩陣。通過聯(lián)合稀疏恢復(fù)理論，支撐集矩陣A可以通過式(19)最優(yōu)化方法獲得

(19)

2.2 字典失配問題

稀疏恢復(fù)STAP方法使用的空時導(dǎo)向矢量字典中空間頻率和多普勒頻率等間隔離散地分布在空時平面上。對于正側(cè)視陣，且折疊系數(shù)β=1時，雜波脊位于空時平面對角線位置，正好位于空時平面網(wǎng)格點上，如圖2(a)所示，字典中空時導(dǎo)向矢量可以準確地與雜波單元對應(yīng)，此時不存在字典失配問題，STAP處理性能較好。當折疊系數(shù)β≠1即雜波脊不在空時平面對角線上，如圖2(b)所示，或工作在非正側(cè)視陣情況，雜波脊在空時平面上呈現(xiàn)曲線形式時，如圖2(c)所示，可以看出雜波單元通常不能與空時導(dǎo)向矢量對應(yīng)，字典失配將嚴重影響雜波譜稀疏恢復(fù)性能。

圖2 雜波脊在空時平面上的分布示意圖

3 字典失配情況下稀疏貝葉斯學習空時譜估計

3.1 空時二維動態(tài)字典模型

(20)

利用空時二維動態(tài)字典估計字典失配誤差Δfd和Δfs后，可得到失配補償后的字典Φ=[Φ1,Φ2,…,ΦNsNd]，即

Φ(Δfd,Δfs)=Ψ+Fdiag(Δfd)+Gdiag(Δfs)

(21)

式中：Ψ為式(12)定義的失配字典。

(22)

(23)

(24)

(25)

估計字典失配誤差進行補償后，式(15)改寫為

x=Φ(Δfd,Δfs)a+n

(26)

式(18)可以改寫為

X=Φ(Δfd,Δfs)A+N

(27)

3.2 字典失配誤差估計及空時譜稀疏恢復(fù)

本節(jié)以式(27)表示的MMV模型，給出字典失配誤差估計及空時譜估計的稀疏貝葉斯學習算法，式(26)表示的SMV模型可由MMV模型方法中快拍數(shù)設(shè)置K=1得到。

噪聲向量N假設(shè)為零均值復(fù)高斯白噪聲，則

(28)

p(X|A,Δfd,Δfs,λ0)=

(29)

式中：λ0未知，假設(shè)服從伽馬先驗分布，即

p(λ0|c0,d0)=Γ(λ0|c0,d0)

(30)

為了保證獲得廣泛的超先驗，一般假設(shè)c0→0，d0→0，例如取c0=1×10-4，d0=1×10-4。

稀疏支撐集A=[a1,a2,…,aK]的各列獨立，且假設(shè)服從復(fù)高斯先驗分布，即

(31)

式中：Λ=diag(λ)，λ假設(shè)為NsNd維獨立伽馬分布，即

(32)

式中：超參數(shù)ρ=1×10-3。

字典失配誤差Δfd和Δfs假設(shè)服從NsNd維均勻分布，即

Δfd～U([-rd/2,rd/2]NsNd)

(33)

Δfs～U([-rs/2,rs/2]NsNd)

(34)

式中：rd和rs分別為離散字典的歸一化多普勒頻率間隔和歸一化空間頻率間隔。

由于后驗分布p(A,Δfd,Δfs,λ0,λ|X)不能顯式給出，將A當作隱變量，則其后驗分布滿足：

(35)

式中:

μk=λ0ΣΦHxk

(36)

Σ=(λ0ΦHΦ+Λ-1)-1

(37)

由式(35)可知，信號支撐集矩陣A的稀疏解與μk和Σ的稀疏解一一對應(yīng)。計算μk和Σ的稀疏解需要對超參數(shù)Δfd，Δfs，λ0和λ進行估計。根據(jù)稀疏貝葉斯理論，超參數(shù)可以用最大后驗方法估計，即最大化p(Δfd,Δfs,λ0,λ|X)，由于p(X)與上述超參數(shù)無關(guān)，因此最大化p(Δfd,Δfs,λ0,λ|X)與最大化p(Δfd,Δfs,λ0,λ,X)=p(Δfd,Δfs,λ0,λ|X)p(X)等價，并且有

p(A,X,Δfd,Δfs,λ0,λ)=

p(X|A,Δfd,Δfs,λ0)p(A|λ)p(λ)p(λ0)·

p(Δfd)p(Δfs)

(38)

利用EM方法迭代更新超參數(shù)λ和λ0，具體推導(dǎo)過程可參見文獻[24]，λ和λ0的更新方程分別為

n=1,2,…,NsNd

(39)

(40)

tr{(Ψ+Fdiag(Δfd)+Gdiag(Δfs))Σ(Ψ+Fdiag(Δfd)+Gdiag(Δfs))H}=

(Δfd)TP1Δfd-2V1Δfd+(Δfs)TP2Δfs-2V2Δfs+C

(41)

式中：C為常數(shù)。

P1=

(42)

P2=

(43)

V1=

(44)

V2=

(45)

最小化式(41)可得到Δfd和Δfs的估計值

2V1Δfd+(Δfs)TP2Δfs-2V2Δfs}

(46)

式(46)分別對Δfd和Δfs的偏導(dǎo)數(shù)為零，則可以得到字典失配誤差的更新，即

(47)

(48)

通過式(21)獲得字典失配誤差補償后的空時導(dǎo)向矢量字典Φ(Δfd,Δfs)，并重復(fù)迭代λ、λ0、Δfd和Δfs的EM估計過程直至達到設(shè)置的估計精度。在兼顧超參數(shù)估計精度和迭代收斂速度的情況下，設(shè)置超參數(shù)誤差預(yù)設(shè)值最大值1×10-3，最大迭代次數(shù)2 000次。

迭代完成后獲得稀疏解A及字典失配誤差，進而獲得修正后的空時導(dǎo)向矢量字典Φ，并重構(gòu)雜波協(xié)方差矩陣，即

(49)

通過式(49)獲得CCM估計值后，雜波空時譜可以表示為

(50)

綜上所述，字典失配情況下的稀疏貝葉斯學習雜波空時譜估計方法具體步驟如下：

步驟1根據(jù)式(12)初始化空時導(dǎo)向矢量字典Ψ。

步驟2根據(jù)式(21)構(gòu)建失配字典動態(tài)模型。

根據(jù)式(36)和式(37)計算均值μk和方差Σ。

根據(jù)式(39)、式(40)、式(47)和式(48)更新超參數(shù)λ、λ0、Δfd和Δfs。

判斷誤差是否達到預(yù)設(shè)值或迭代次數(shù)是否達到最大值，不滿足則重復(fù)步驟2。

步驟3根據(jù)式(49)估計雜波協(xié)方差矩陣，然后根據(jù)式(50)計算雜波空時譜。

4 仿真實驗

為了驗證本文方法的有效性，通過如下參數(shù)生成仿真數(shù)據(jù)。雷達工作頻率fo=450 MHz，載機平臺高度H=9 000 m，陣元數(shù)M=8，相干脈沖數(shù)N=8，脈沖重復(fù)頻率fr=300 Hz，360個雜波單元在0°～180°均勻分布，雜噪比CNR=40 dB， 從距離R0=15 km處開始仿真3個距離單元的雜波數(shù)據(jù)，距離分辨率為37.5 m。實驗對比了本文方法和使用固定離散字典稀疏貝葉斯學習的SBL-STAP方法[18]。SBL-STAP方法中最大迭代次數(shù)為1 000次，正則化參數(shù)初始值取1×10-2，超參數(shù)修剪閾值取1×10-4。

4.1 雜波空時功率譜恢復(fù)效果比較

實驗1為了比較不同偏航角情況(有無字典失配問題)對稀疏恢復(fù)STAP的影響，陣列采用正側(cè)視模式。折疊系數(shù)β=1，即雜波脊的斜率為1，對應(yīng)載機運動速度vp=50 m/s。空時導(dǎo)向矢量字典網(wǎng)格劃分系數(shù)ρs=ρd=4，即歸一化空間頻率和歸一化多普勒頻率都被等間隔分為32份，這時雜波脊正好落在空時導(dǎo)向矢量字典網(wǎng)格點上，不存在字典失配的問題。圖3為真實雜波譜及上述兩種方法使用相同的3個雜波樣本進行稀疏恢復(fù)得到的雜波譜，可以看出當不存在字典失配時，SBL-STAP方法和本文方法形成的雜波譜均集中在雜波脊上，無明顯展寬，具有較好的雜波抑制性能。

圖3 不存在字典失配問題時的雜波譜估計

實驗2陣列采用正側(cè)視模式，折疊系數(shù)為β=0.8，即雜波脊的斜率為0.8，對應(yīng)載機運動速度vp=40 m/s?？諘r導(dǎo)向矢量字典網(wǎng)格仍為實驗1中ρs=ρd=4的等間隔劃分，由于雜波脊斜率與字典網(wǎng)格不匹配，雜波脊不能完全落在空時導(dǎo)向矢量字典網(wǎng)格上，此時存在字典失配問題，如圖4(a)所示。從圖4(b)中可以看出， SBL-STAP方法恢復(fù)的雜波譜有明顯的展寬，分辨率不高，雜波抑制性能損失較大。由實驗1和實驗2的結(jié)果對比可以看出，這種雜波譜展寬的情況主要是由空時導(dǎo)向矢量字典失配造成的，說明字典網(wǎng)格失配會使得稀疏恢復(fù)性能明顯下降，STAP處理效果較差。而本文提出的基于動態(tài)字典的稀疏貝葉斯學習空時自適應(yīng)處理方法，利用動態(tài)字典實現(xiàn)字典失配誤差的估計補償，能夠形成較為清晰的雜波脊、連續(xù)性和展寬情況都獲得了比較理想的結(jié)果，如圖4(c)所示。

圖4 存在字典失配問題時的雜波譜估計(正側(cè)視陣，β=0.8)

實驗3陣列采用非正側(cè)視陣，φ=30°，載機運動速度vp=50 m/s真實雜波如圖5(a)所示。從圖5(b)可以看出SBL-STAP方法恢復(fù)得到的雜波譜中雜波脊較實驗2中情況更為模糊，這主要是由于非正側(cè)視陣的雜波脊為橢圓形，存在更多不能落在空時平面網(wǎng)格點上的雜波單元，由于字典失配問題的影響，稀疏恢復(fù)得到的雜波譜誤差較大。而如圖5(c)所示，本文方法獲得的雜波譜仍然能夠形成較為清晰的雜波脊，在非正側(cè)視陣存在字典失配的情況下較SBL-STAP方法有較大的優(yōu)勢。

4.2 雜波抑制性能比較

除上述雜波譜恢復(fù)情況分析外，STAP處理的雜波抑制效果通常使用信雜噪比損失(Signal to Clutter plus Noise Ratio Loss, SCNRLoss)來衡量[26]，即

(51)

式中：

(52)

圖6為進行100次Monte Carlo實驗得到的平均結(jié)果。其中圖6(a)為正側(cè)視陣情況，折疊系數(shù)β=0.8存在字典失配情況下的信雜噪比損失。圖中OPT對應(yīng)曲線表示利用CCM真實值進行STAP處理獲得的信雜噪比損失，稱為最優(yōu)情況?？梢钥闯霰疚姆椒ㄐ纬闪溯^窄的零陷，在抑制雜波的同時，能夠獲得較好的運動目標檢測性能，除了在主雜波區(qū)形成零陷外，其他區(qū)域的SCNRLoss均較小，僅比最優(yōu)情況低2 dB左右。而SBL-STAP方法形成的零陷較寬，這會使得更多與雜波多普勒頻率接近的運動目標被當作雜波抑制掉，大大降低系統(tǒng)的運動目標檢測性能。圖6(b)為非正側(cè)視陣，φ=30°存在字典失配的情況下的信雜噪比損失?？梢钥吹脚c圖6(a)相似的結(jié)果，本文方法主雜波區(qū)形成較窄的零陷，在抑制雜波的同時，能夠獲得較好的運動目標檢測性能。而使用固定離散字典的SBL-STAP方法形成的雜波抑制區(qū)域較寬，對于運動目標的檢測性能下降較為嚴重。

圖6 字典失配情況下的信雜噪比損失

5 結(jié) 論

稀疏恢復(fù)STAP方法受字典失配影響較大，雜波抑制性能下降明顯。本文提出了一種利用空時二維動態(tài)字典對失配字典誤差進行估計補償?shù)南∈柝惾~斯學習空時譜估計算法，實現(xiàn)了存在字典失配情況下的雜波譜有效恢復(fù)，能夠有效提高機載雷達雜波抑制性能。仿真實驗表明，本文方法在存在字典失配的情況下，可以有效恢復(fù)雜波協(xié)方差矩陣，提高雜波抑制效果，其性能優(yōu)于固定離散字典的SBL-STAP方法。