郭藝奪，宮 健

(空軍工程大學(xué)防空反導(dǎo)學(xué)院，陜西西安 710051)

0 引言

空時自適應(yīng)處理(Space Time Adaptive Processing, STAP)[1-2]方法能夠抑制地面雜波，檢測慢速目標，其雜波抑制性能取決于由訓(xùn)練樣本單元估計的雜波協(xié)方差的準確度。在實際中，待測單元的雜波協(xié)方差是由與待測單元獨立同分布(Independent Identically Distributed, IID)的訓(xùn)練樣本單元估計得到的，但在機載非正側(cè)陣雷達系統(tǒng)中，雜波多普勒頻率和空間頻率具有嚴重的距離相關(guān)性，尤其是近程條件下，距離相關(guān)性更加顯著，使得訓(xùn)練樣本單元與待測單元波分布不滿足IID條件[3-4]，進而無法以訓(xùn)練樣本單元來準確估計待測單元的雜波協(xié)方差矩陣，使得STAP方法的雜波抑制性能下降，無法有效檢測慢速目標。

補償距離相關(guān)性的方法有很多，包括多普勒彎曲(Doppler Warping, DW)[5]、角度多普勒補償(Angle Doppler Compensation, ADC)[6]、空時內(nèi)插補償(Space Time Interpolating Technique, STINT)[7]和基于配準補償(Registration Based Compensation, RBC)[8]等方法。自適應(yīng)角度多普勒補償(Adaptive Angle Doppler Compensation, A2DC)方法[9]及其改進方法[10]是一種補償非正側(cè)陣雷達雜波距離相關(guān)性的常用方法。該方法利用雷達雜波回波數(shù)據(jù)通過子孔徑平滑估計雜波主特征向量的最小方差無失真響應(yīng)(Minimum Variance Distortionless Response, MVDR)譜，基于此估計訓(xùn)練樣本單元和待測單元的雜波譜中心位置，然后將訓(xùn)練樣本單元的雜波譜中心分別沿角度和多普勒方向平移，使得平移后訓(xùn)練樣本單元的MVDR譜中心與待檢測單元的MVDR譜中心重合。由于該方法利用雷達雜波回波數(shù)據(jù)自適應(yīng)計算各訓(xùn)練樣本單元的轉(zhuǎn)換矩陣，因此可有效降低了載機平臺提供的飛行配置參數(shù)誤差對補償性能造成的影響。

然而，A2DC方法估計的雜波分布特性是基于子孔徑平滑獲得的，孔徑損失使得空時譜分辨率有一定程度的損失，系統(tǒng)自由度(Degree of Freedom, DOF)下降，影響后續(xù)補償過程性能，導(dǎo)致變換后的各訓(xùn)練樣本單元平穩(wěn)性下降，雜波抑制性能受限。此外，A2DC方法在計算各距離單元雜波的主特征向量的MVDR譜中心位置時，通常先要利用特征分解求得各距離單元雜波協(xié)方差矩陣的主特征向量，而典型的矩陣特征分解方法如特征值分解和奇異值分解算法的運算量都非常大，這使得A2DC方法實現(xiàn)起來較為復(fù)雜，限制了其應(yīng)用。文獻[10]利用投影逼近子空間追蹤(Projection Approximation Subspace Tracking，PAST)方法解決A2DC方法運算量大的問題，但PAST方法并不是全局收斂的，可能陷入局部極小點而無法得到全局最優(yōu)解。

針對上述問題，本文提出了基于稀疏恢復(fù)(Sparse Recovery, SR)[11]的A2DC方法——SR-A2DC方法。該方法首先對不同距離單元的回波數(shù)據(jù)進行稀疏恢復(fù)，得到雜波空時分布譜，并計算雜波協(xié)方差矩陣，已避免常規(guī)A2DC方法孔徑損失帶來的影響；然后利用正交近似投影子空間追蹤(Orthogonal PAST, OPAST)[12]方法估計不同距離單元的雜波主特征向量，降低運算量；最后估計主特征向量的空間頻率和多普勒頻率，計算空時轉(zhuǎn)換矩陣，將訓(xùn)練樣本單元回波數(shù)據(jù)與待測單元的回波數(shù)據(jù)進行配準。理論分析和仿真實驗表明，該方法可以實現(xiàn)雜波距離相關(guān)性的自適應(yīng)補償，提高主瓣雜波的補償效果，減少運算量，在存在誤差時也能獲得較好的處理性能。

1 雜波模型及A2DC方法

1.1 雜波模型

機載雷達第l個距離單元的雜波數(shù)據(jù)由該距離單元上的多個離散雜波塊疊加而成：

(1)

式中，P為離散雜波塊個數(shù)，σi為第i個雜波塊對應(yīng)的復(fù)幅度，wt和ws分別為雜波的歸一化多普勒頻率和空間頻率，Si為對應(yīng)的空時二維導(dǎo)向矢量：

Si(wt,ws)=St,i(wt,i)?Ss,i(ws,i)

(2)

式中，St,i和Ss,i分別為時域?qū)蚴噶亢涂沼驅(qū)蚴噶浚?/p>

(3)

式中，wt,i=2πfd,i/fprf，ws,i=2πdfs,i/λ，N和K分別為天線陣元數(shù)和相干脈沖數(shù)，fd,i和fs,i分別為雜波對應(yīng)的多普勒頻率和空間頻率，d和λ分別為陣元間距和波長。

非正側(cè)陣雷達雜波的多普勒頻率和空間頻率與俯仰角和方位角的耦合關(guān)系為

(4)

fs,i=cosφicosθi

(5)

式中，θi和φi分別為雜波對應(yīng)的方位角和俯仰角，v0為載機速度，θp為天線陣面與載機速度之間的夾角。

1.2 A2DC方法

A2DC方法首先以空域和時域的孔徑損失為代價估計各距離單元的雜波協(xié)方差矩陣；接著對估計得到的雜波協(xié)方差矩陣進行特征分解，得到主特征向量，然后計算各個距離單元主特征向量的MVDR譜中心位置；最后以待測單元為基準，通過角度和多普勒平移，使得各訓(xùn)練樣本單元主特征向量的MVDR譜中心與待測單元主特征向量的MVDR譜中心重合。其具體實現(xiàn)步驟如下。

第一步，利用子孔徑平滑法獲得該距離單元的雜波協(xié)方差矩陣。

子孔徑平滑法首先將第l個距離單元的雷達回波數(shù)據(jù)Xl∈CNK×1變換成矩陣Xl∈CN×K：

(6)

設(shè)空域子孔徑和時域子孔徑分別為P和Q，對Xl進行空時平滑可以得到D=(N-P+1)·(K-Q+1)個子矩陣Xu,v∈CP×Q為

(7)

式中：u=1,2,…,N+P-1；v=1,2,…,K-Q+1。

由子孔徑平滑得到的雜波回波數(shù)據(jù)子矩陣，可估計第l個距離單元回波的協(xié)方差矩陣Rsl∈CPQ×PQ為

(8)

其中，Vec(Xu,v)表示將矩陣Xu,v的第2列置于第1列的下面，第3列置于第2列的下面，以此類推，使矩陣Xu,v變?yōu)橐涣邢蛄?。D的值應(yīng)不小于雜波協(xié)方差矩陣Rl的列數(shù)，即D≥PQ，以確保Rl為滿秩矩陣。

第二步，對雜波協(xié)方差矩陣進行特征值分解，得到相應(yīng)的主特征向量。

在得到各距離單元的雜波協(xié)方差矩陣估計后，A2DC方法利用特征值分解求得雜波協(xié)方差矩陣Rl的最大特征值λl,max及其對應(yīng)的最大特征向量hl,max。

第三步，根據(jù)主特征向量的相位信息估計得到主特征向量的MVDR譜中心的多普勒頻率和空間頻率。

主特征向量的MVDR譜中心的位置可以通過二維峰值搜索計算得到，也可以直接通過主特征向量的相位信息估計得到。后者需要估計雜波協(xié)方差矩陣Rl的主特征向量的空間相位斜率κl,s和時間相位斜率κl,t，利用它們與空間頻率和多普勒頻率的線性關(guān)系計算主特征向量的MVDR譜中心的多普勒頻率fl,d和空間頻率fl,s。

(9)

第四步，計算訓(xùn)練單元與待測單元主特征向量MVDR譜中心的多普勒頻率和空間頻率之間的差值為

(10)

接著，根據(jù)訓(xùn)練樣本單元和待測單元的MVDR譜中心的多普勒頻率和空間頻率，計算各個訓(xùn)練樣本單元的空時轉(zhuǎn)換矩陣TA2DC,l，對訓(xùn)練樣本單元接收回波數(shù)據(jù)進行補償，其中，第l個訓(xùn)練單元的轉(zhuǎn)換矩陣為

TA2DC,l=Ttl?Tsl

(11)

式中，

(12)

(13)

經(jīng)過A2DC補償后的第l個距離單元的雷達回波數(shù)據(jù)為

(14)

2 SR-A2DC方法

由以上分析可知，A2DC方法可以自適應(yīng)實現(xiàn)雜波距離相關(guān)性補償，但該方法需要對雜波協(xié)方差矩陣Rl進行特征分解，使得計算空時轉(zhuǎn)換矩陣TA2DC,l的運算量巨大，且由于孔徑損失，Rl的估計準確度不高，以上兩點使得A2DC方法的應(yīng)用受限。

2.1 基于稀疏的雜波協(xié)方差矩陣求解

由式(1)可知，機載雷達回波數(shù)據(jù)是由不同空間頻率和多普勒頻率的回波數(shù)據(jù)疊加而成的，將空間頻率和多普勒頻率分別遍歷并離散為Ns=ρsN，Nd=ρdN個分辨單元，則第l個距離單元的雜波回波數(shù)據(jù)可以表示為

(15)

式中，ρs和ρd分別表示空間頻率和多普勒頻率的離散化程度，在高分辨情況下一般遠大于1；Sq(wt,q,ws,q)為第q個空時導(dǎo)向矢量，γq為其對應(yīng)的復(fù)幅度，wt,q和ws,q為對應(yīng)的多普勒頻率和空間頻率；αl=[αl,1,αl,2,…,αl,NsNd]代表雷達回波數(shù)據(jù)在空間頻率-多普勒頻率域上的幅度分布，即雜波空時譜；Ψ為超完備基矩陣：

(16)

估計雜波空時譜等同于在方程(15)中已知Xk和Ψ而求解αl。由于Ψ的列數(shù)NsNd遠大于行數(shù)NK，因此方程(15)屬于欠定方程，存在多個可能解。根據(jù)稀疏恢復(fù)理論，當雜波空時譜αl具有稀疏性時，欠定方程可以較高的概率求解。

針對上述問題的求解方法主要有3種，即凸優(yōu)化方法[13]、FOCUSS算法[14]和MP算法[15]，其中，F(xiàn)OCUSS算法不利用可能存在誤差的先驗知識，采用后驗知識迭代加權(quán)，使解的能量逐步集中，可以避免誤差帶來的不利影響，且FOCUSS算法在初值設(shè)置合理時通?？色@取較好的譜估計性能，同時運算量相比范數(shù)最小化方法大大降低。因此綜合考慮算法性能和運算量，本文選擇FOCUSS算法進行計算。

實際中，僅利用待測單元恢復(fù)得到的雜波空時譜與真實雜波空時譜相比存在一定差異，一方面，如果待測單元中存在目標，則后續(xù)的處理會將目標同時抑制；另一方面，稀疏恢復(fù)估計的雜波空時譜存在過于稀疏的問題。為了更準確地表示雜波空時譜分布，一般采用對臨近訓(xùn)練樣本單元稀疏恢復(fù)進行平均的方法。需要特別強調(diào)的是，訓(xùn)練樣本單元必須選擇與待測單元相鄰，以保證雜波近似獨立同分布。

對式(15)進行求解，得到雜波的空時譜分布αl，則可獲得對應(yīng)的雜波協(xié)方差矩陣估計：

Sq(wt,q,ws,q)H

(17)

式中，αq為雜波空時譜估計在第q個位置上的雜波復(fù)幅度；Sq(wt,q,ws,q)為第q個空時導(dǎo)向矢量。

2.2 引入OPAST算法估計雜波主特征向量

A2DC方法的第二步是對雜波協(xié)方差矩陣進行分解得到主特征向量，但是傳統(tǒng)的特征分解方法如特征值分解和奇異值分解等，雖然性能優(yōu)越，但運算量大，特別是對于大線陣而言，不利于實時實現(xiàn)?？焖僮涌臻g跟蹤方法是一種新的特征值分解方法，其中PAST和PASTd算法收斂速度快，誤差小，但該方法不是全局收斂的，而正交PAST(OPAST)算法相比PAST算法，在運算量增加不多的情況下，可以保證全局收斂。本文利用OPAST算法對雜波協(xié)方差矩陣進行特征值分解，獲得主特征向量，以降低A2DC方法的運算量。

OPAST算法將特征子空間的確定轉(zhuǎn)化為一種無約束的最優(yōu)化問題進行求解，定義的無約束代價函數(shù)為

(18)

式中，Rl為第l個距離單元的雜波協(xié)方差矩陣估計，矩陣變元Wl∈NK×r(r≤NK)，Tr(·)表示矩陣的跡。

在利用遞推最小二乘算法求解代價函數(shù)的全局最小值時，為保證全局收斂和矩陣變元Wl的規(guī)范正交性，在每次迭代中對矩陣變元Wl進行正交化，即引入正交化公式：

(19)

當代價函數(shù)J(Wl)取得最小值時，Wl收斂于雜波子空間的一組正交基上，特別地，當r=1，即Wl∈NK×1時，最小化J(Wl)得到Wl，即為Rl的歸一化最主要特征向量。

與常規(guī)A2DC方法相比，利用OPAST算法求解雜波協(xié)方差矩陣的主特征向量，運算量從O(NK3)減少為4NKS+O(S2)，其中，S為雜波協(xié)方差矩陣Rl的秩。

3 仿真實驗

選取斜側(cè)陣(θp=60°)機載相控陣雷達進行仿真實驗，其中實驗條件為：發(fā)射接收陣元數(shù)均為12；相干脈沖間隔內(nèi)的脈沖數(shù)為12；陣元間隔0.16 m，載機高度8 km，載機速度120 m/s，波長0.32 m，脈沖重復(fù)頻率3 000 Hz；雜噪比60 dB，雷達最大作用距離800 km；雷達主波束方位角為90°；待測單元145，假設(shè)對應(yīng)的距離為25 km，且不存在距離模糊。

采樣協(xié)方差求逆方法(LSMI)的訓(xùn)練樣本數(shù)為2NK；稀疏恢復(fù)中，訓(xùn)練樣本數(shù)為4；空間頻率和多普勒頻率離散化程度均為8，即ρs=ρd=8；子空間平滑的空域子孔徑和時域子孔徑為P=Q=6。

實驗1：雜波空時譜估計

本實驗對比真實雜波功率譜、未經(jīng)補償?shù)腖SMI法的雜波功率譜、常規(guī)A2DC方法子孔徑平滑法的雜波空時譜和本文SR-A2DC方法估計的雜波空時譜，如圖1～圖4所示。

圖1 真實雜波功率譜

圖2 未經(jīng)補償?shù)腖SMI法估計的雜波功率譜

圖3 A2DC方法子孔徑雜波空時譜

圖4 SR-A2DC方法全孔徑雜波空時譜

由圖1和圖2可以看出，在非正側(cè)陣的情況下，由于雜波具有距離相關(guān)性，使用LSMI方法直接利用訓(xùn)練樣本單元估計待測單元的雜波協(xié)方差矩陣，不同距離單元的雜波疊加在一塊，會導(dǎo)致雜波譜嚴重展寬。由圖3和圖4可以看出，常規(guī)A2DC方法得到的雜波空時譜，由于子孔徑平滑造成的孔徑損失，存在能量擴散和部分雜波分布難以估計的問題；SR-A2DC方法利用FOCUSS方法估計雜波空時譜，不存在孔徑損失，因此準確度和分辨率更高。

實驗2：雜波譜中心多普勒頻率和空間頻率

本實驗對比常規(guī)A2DC方法特征值分解(Eigen Decomposition, ED)和本文SR-A2DC方法OPAST算法估計的雜波譜中心多普勒頻率和空間頻率與真實值之間的差異，如圖5、圖6所示。

圖5 歸一化多普勒頻率對比

圖6 歸一化空間頻率對比

從圖5和圖6可以看出，相比常規(guī)A2DC方法特征值分解估計的雜波譜中心多普勒頻率和空間頻率，本文SR-A2DC方法OPAST算法準確性更高。這是因為，OPAST算法能夠以較高的收斂速度保證較小的誤差，而基于稀疏恢復(fù)的雜波空時譜估計不存在孔徑損失，兩者結(jié)合性能優(yōu)于常規(guī)A2DC方法。

實驗3：補償效果對比

本實驗對比常規(guī)A2DC方法和SR-A2DC方法對雜波距離相關(guān)性進行補償?shù)难a償效果，并以改善因子為基準衡量不同方法的雜波抑制性能，如圖7～圖9所示。其中，改善因子的定義為輸出信雜噪比(Signal to Clutter-plus-Noise Ratio, SCNR)與輸入SCNR的比值。

由圖7和圖8可知，A2DC方法和SR-A2DC方法均能對非正側(cè)陣雜波距離相關(guān)性進行補償，經(jīng)過補償后，主瓣雜波基本配準，但由于A2DC方法是基于子孔徑平滑，存在孔徑損失，補償性能在一定程度上受限，而SR-A2DC方法是基于全孔徑的，不存在孔徑損失，可獲得更高的分辨率，設(shè)計的空時轉(zhuǎn)換矩陣更加有效，因此補償性能更優(yōu)。

圖7 A2DC方法雜波距離相關(guān)性補償效果

圖8 SR-A2DC方法雜波距離相關(guān)性補償效果

由圖9可知，經(jīng)過本文SR-A2DC方法進行雜波距離相關(guān)性補償后，利用LSMI方法估計雜波協(xié)方差矩陣，計算空時濾波器權(quán)值，可以在主雜波處形成深凹口，且相比A2DC方法更窄，對慢速目標的檢測能力可提高5～10 dB。而利用A2DC方法設(shè)計的空時濾波器，則存在著凹口偏移、展寬等問題。

圖9 雜波抑制性能

4 結(jié)束語

本文詳細分析了常規(guī)A2DC方法的原理和步驟，針對該方法存在孔徑損失和運算量大的問題，提出了基于稀疏恢復(fù)的SR-A2DC方法，并將OPAST算法引入到雜波協(xié)方差矩陣特征分解中。理論分析和仿真實驗表明：1)所提的SR-A2DC方法可以避免空時損失帶來的影響，減少運算量；2)SR-A2DC方法與常規(guī)A2DC方法，雜波距離相關(guān)性補償性能更優(yōu)，更適于工程實現(xiàn)；3)由于本文使用稀疏恢復(fù)估計雜波空時譜，而稀疏恢復(fù)求解所構(gòu)造的超完備矩陣存在著基失配(off-grid)等問題，因此，如何減少或者消除off-grid的影響，自適應(yīng)設(shè)計超完備基矩陣，獲得更加準確的雜波空時譜分布，是下一步研究工作的重點。