一種基于聯(lián)合稀疏恢復(fù)的空時自適應(yīng)處理方法

段克清 王澤濤 謝文沖 高 飛 王永良

①(空軍預(yù)警學(xué)院 武漢 430019)

②(國防科學(xué)技術(shù)大學(xué)電子科學(xué)與工程學(xué)院 長沙 410073)

1 引言

機(jī)載相控陣?yán)走_(dá)體制下，空時自適應(yīng)處理(Space-Time Adaptive Processing, STAP)可實(shí)現(xiàn)對強(qiáng)雜波的有效抑制以及慢速運(yùn)動目標(biāo)的檢測[1,2]。該技術(shù)需利用相鄰多個距離門數(shù)據(jù)(訓(xùn)練樣本)估計(jì)待檢測距離門雜波協(xié)方差矩陣(Clutter Covariance Matrix, CCM)，進(jìn)而計(jì)算空時自適應(yīng)濾波器權(quán)系數(shù)。根據(jù)RMB (Reed-Mallett-Brennan)準(zhǔn)則，為使得由CCM估計(jì)誤差所帶來的雜波抑制性能損失小于3 dB，所需獨(dú)立同分布訓(xùn)練樣本數(shù)至少為2倍系統(tǒng)自由度[3]。在實(shí)際應(yīng)用中，由于機(jī)載雷達(dá)往往面臨復(fù)雜多變的下視場景，常常難以獲取足夠的均勻訓(xùn)練樣本，因此小樣本、高性能STAP算法是目前該領(lǐng)域中的熱點(diǎn)研究方向。

近年來，壓縮感知或稀疏恢復(fù)技術(shù)被廣泛應(yīng)用于信號處理領(lǐng)域，并初步應(yīng)用于機(jī)載雷達(dá)STAP技術(shù)中[4－10]。由于天線的方向性，雜波在空時 2維譜上的分布是稀疏的，因此可利用少量甚至單個訓(xùn)練樣本基于超完備字典來獲取空時譜上強(qiáng)雜波點(diǎn)的位置，即稀疏系數(shù)中的非零元素位置。超完備字典可利用超分辨空時2維導(dǎo)向矢量構(gòu)造，求解稀疏系數(shù)的問題實(shí)則為稀疏約束下欠定方程求解問題。針對該問題的求解目前主要有3種途徑，即凸優(yōu)化方法[11－13]、貪婪算法[14－16]和 FOCUSS 算法[17]。其中，貪婪算法運(yùn)算量較小，算法性能比較穩(wěn)定，應(yīng)用較為廣泛，因此本文重點(diǎn)討論貪婪算法在STAP技術(shù)中的應(yīng)用。目前，文獻(xiàn)[6-8]主要探討了利用單樣本恢復(fù)或多個樣本分別恢復(fù)后再聯(lián)合處理的稀疏恢復(fù)STAP算法，其本質(zhì)上仍為單觀測矢量情況(Single Measurement Vector, SMV)；文獻(xiàn)[9,10]針對稀疏恢復(fù)STAP問題進(jìn)行了深入且系統(tǒng)的研究，一方面研究了基于稀疏濾波器的STAP技術(shù)，另一方面針對空時譜稀疏問題重點(diǎn)研究了非參數(shù)化迭代稀疏恢復(fù)算法，旨在增強(qiáng)算法對正則化參數(shù)設(shè)置的穩(wěn)健性和運(yùn)算量的降低。然而，上述研究均未深入涉及多觀測矢量(Multiple Measurement Vector, MMV)聯(lián)合稀疏問題。

實(shí)際上，近年來在腦成像等領(lǐng)域的技術(shù)需求，激發(fā)了人們基于MMV進(jìn)行聯(lián)合稀疏恢復(fù)的興趣[18]，也提出了一些較為有效的算法[18－23]。其中同時正交匹配追蹤算法(Simutaneous Orthogonal Matching Pursuit, SOMP)[18－21]是貪婪算法由SMV到MMV的直接拓展，具有較低的運(yùn)算量，同時在無噪聲背景下具有理想的信號恢復(fù)性能，是目前最為高效和典型的聯(lián)合稀疏恢復(fù)算法。然而，在實(shí)際應(yīng)用背景下，機(jī)載雷達(dá)回波信號中除目標(biāo)和強(qiáng)雜波外不可避免地存在各種噪聲，因此，如果簡單將SOMP算法引入到STAP技術(shù)中，受噪聲影響，雜波抑制性能的穩(wěn)定性將得不到保證。

針對上述問題，本文提出了一種基于雜波子空間的貪婪算法進(jìn)行聯(lián)合稀疏恢復(fù)，并應(yīng)用于STAP進(jìn)行雜波抑制處理。首先，該方法充分利用了多個訓(xùn)練樣本中的雜波信息，因此較現(xiàn)有稀疏恢復(fù)STAP方法具備更強(qiáng)的信號恢復(fù)能力；其次，不同于 SOMP算法中利用剩余雜波矩陣來尋找最大相關(guān)原子的策略，該方法采用了剩余雜波子空間投影來尋找最大相關(guān)原子，并利用2l范數(shù)作為準(zhǔn)則來衡量相關(guān)性的大小，從而有效回避了雜波矩陣中的噪聲分量對稀疏恢復(fù)算法的不利影響，同時可高效地尋找到過完備字典中相關(guān)性最強(qiáng)的原子。本文最后利用仿真實(shí)驗(yàn)驗(yàn)證了所提STAP算法的有效性，并給出了相應(yīng)結(jié)論。

2 機(jī)載雷達(dá)信號模型及空時譜估計(jì)

不考慮距離模糊情況下，機(jī)載雷達(dá)在某個距離門上的回波信號主要由該距離環(huán)上多個離散雜波塊的回波疊加而成：

其中p為距離環(huán)上雜波塊個數(shù)，γp為第p個雜波塊的散射系數(shù)，s為目標(biāo)信號，n為噪聲，?p為第p個雜波塊對應(yīng)空時2維導(dǎo)向矢量：

其中, St,p和 Ss,p分別為第p個雜波塊對應(yīng)的時間導(dǎo)向列矢量和空間導(dǎo)向列矢量，且有

其中，θs,p為第p個雜波塊對應(yīng)的空間錐角，fd,p=2V cos( θs,p)/λ 為第p個雜波塊對應(yīng)多普勒頻率，V為載機(jī)速度，λ為波長，N為天線陣元數(shù)，d為陣元間距，M為相干脈沖數(shù)，fr為脈沖重復(fù)頻率。

由式(1)可知，機(jī)載雷達(dá)回波信號實(shí)際上由不同來向不同多普勒頻率的回波信號疊加而成，如果將所有空間頻率和多普勒頻率均遍歷并分別離散化為Ns=ρsN和Nd=ρdM個頻率點(diǎn)，且只考慮訓(xùn)練樣本中回波信號，則式(1)中信號還可表示為：

其中，ρs和ρd分別表示離散化程度，在高分辨情況下一般均遠(yuǎn)大于1, αq為第q個空時網(wǎng)格單元幅度，?q為第q個空時網(wǎng)格單元對應(yīng)的空時2維導(dǎo)向矢量，。由于字典ψ的行數(shù)NM遠(yuǎn)小于列數(shù) NsNd，因此式(5)為病態(tài)方程或欠定方程。

即從式(6)中求解稀疏系數(shù)矩陣A中的非零元素位置。需要注意的是，本文研究背景為正側(cè)視機(jī)載雷達(dá)，因此各距離門中雜波分布特性一致，也就是說各距離門中雜波的稀疏結(jié)構(gòu)可認(rèn)為是相同的。式(6)最終需恢復(fù)得到稀疏系數(shù)矩陣A中非零行的坐標(biāo)，即支撐集。

一旦求得A，則可計(jì)算出角-多普勒網(wǎng)格上雜波所在位置網(wǎng)格點(diǎn)的幅度，得到雜波空時2維譜的高分辨估計(jì)，進(jìn)而可構(gòu)造空時自適應(yīng)權(quán)值進(jìn)行雜波抑制處理。

3 基于子空間正交匹配追蹤的STAP方法(Subspace-based OMP, SbOMP)

3.1 算法原理

無噪聲情況下，貪婪算法可有效地尋找到支撐集。然而，當(dāng)觀測矢量中含有噪聲分量時，現(xiàn)有貪婪算法的稀疏恢復(fù)性能會出現(xiàn)較大波動。針對該問題，本文提出了一種基于子空間的貪婪算法進(jìn)行聯(lián)合稀疏恢復(fù)，并利用其估計(jì)空時2維譜構(gòu)造自適應(yīng)權(quán)值進(jìn)行雜波抑制處理。該方法對剩余雜波矩陣進(jìn)行SVD或QR分解，得到剩余雜波子空間，再構(gòu)造投影矩陣基于l2范數(shù)匹配尋找字典中最大相關(guān)原子，然后將該原子對應(yīng)觀測矩陣(訓(xùn)練樣本)中雜波分量剔除，得到新的剩余雜波矩陣，再返回第1步，不斷迭代，直到滿足終止條件。在通過聯(lián)合稀疏恢復(fù)得到空時譜中強(qiáng)雜波分量位置后，基于最小二乘可計(jì)算出空時譜中雜波分量幅度，進(jìn)而估計(jì)得到雜波協(xié)方差矩陣，最終形成空時2維自適應(yīng)權(quán)值進(jìn)行雜波抑制處理。該方法具體表述如下：

輸入：

?NM× L 維訓(xùn)練樣本矩陣X。

?稀疏度K。

輸出：

?含有K個坐標(biāo)的集合ΛK。

?NM× NM維投影矩陣 PK。

?NM× L維近似矩陣 ZK。

?NM× L維剩余雜波矩陣 YK。

?NM× NM維雜波協(xié)方差矩陣R。

? NM× 1維空時導(dǎo)向矢量W。

程序：

第 1步：初始化剩余雜波矩陣Y0=X，坐標(biāo)集Λ0=?和迭代計(jì)數(shù)k=1。

第2步：計(jì)算第k次迭代后剩余雜波矩陣對應(yīng)的雜波子空間Sk?1以及該子空間對應(yīng)的投影算子Pk?1。

第3步：通過求解以下最優(yōu)化問題求解第k次迭代時的坐標(biāo)kλ

第4步：將第k次迭代時的坐標(biāo)λk并入坐標(biāo)集Λ

第5步：根據(jù)當(dāng)前坐標(biāo)集Λk對應(yīng)字典中的原子計(jì)算投影算子 Pk。

第6步：獲取新的雜波分量 Zk并根據(jù)該分量求得剩余雜波分量 Yk

第7步：令 k = k+ 1并返回第2步直到k = K。

第8步：計(jì)算雜波協(xié)方差矩陣R和空時自適應(yīng)權(quán)值W

第2步中可通過SVD分解或QR分解算法[24]來計(jì)算雜波子空間；第3步中，? Λk表示坐標(biāo)集?中除Λk內(nèi)所含坐標(biāo)以外的所有坐標(biāo)的集合；第7步中，迭代次數(shù)K與雜波秩有關(guān)，根據(jù)經(jīng)驗(yàn)一般選取K為1.5～2倍雜波秩r；第8步中，ψΛ為由ψ中相對于Λ的原子組成的矩陣，AΛ為相對于Λ的稀疏系數(shù)矩陣，σ2為噪聲方差，目標(biāo)空時導(dǎo)向矢量ψ0可通過式(2)求取。

3.2 方法分析

已有貪婪算法均利用剩余雜波矩陣與各原子進(jìn)行相關(guān)，依據(jù)相關(guān)性大小來篩選出與雜波分量最大相關(guān)的原子。無噪聲情況下，該類方法可高效篩選出r個與雜波相關(guān)性最強(qiáng)的原子。然而，存在噪聲時受噪聲分量影響，該類方法每次迭代并非恰好篩選出與雜波相關(guān)的原子，如果迭代次數(shù)依舊為r個，則篩選出與雜波相關(guān)的原子往往是欠缺的，因此一般迭代次數(shù)選取比r稍大來確保所有與雜波相關(guān)的原子均被搜索到。同時，與噪聲分量相關(guān)的部分原子不可避免地也被搜索到，表現(xiàn)為雜波譜上出現(xiàn)部分偽峰，導(dǎo)致雜波抑制性能下降。

本文方法利用剩余雜波子空間構(gòu)造投影算子來篩選與雜波分量相關(guān)的原子，保證了每次與原子進(jìn)行相關(guān)運(yùn)算的分量為純雜波分量，因此不會受到噪聲等其他因素的干擾。采用該算法每次篩選到的原子必然分布在理想雜波脊上，并隨著迭代次數(shù)的增加連續(xù)性增強(qiáng)。如果迭代次數(shù)太少，則體現(xiàn)在雜波譜上較為離散；而如果迭代次數(shù)太多，一方面運(yùn)算量較大；另一方面會導(dǎo)致自動搜索到雜波脊附近相關(guān)性相對較強(qiáng)的原子；因此，根據(jù)經(jīng)驗(yàn)值，一般選取迭代次數(shù)為1.5～2倍雜波秩r。

4 仿真實(shí)驗(yàn)

選取正側(cè)視均勻線陣機(jī)載相控陣?yán)走_(dá)進(jìn)行仿真實(shí)驗(yàn)，其中載機(jī)高度為 9000 m，載機(jī)速度為 150 m/s，波長為0.3 m，重頻為2000 Hz，主波束俯仰角為3o，天線陣元個數(shù)N=8，相干脈沖數(shù)M=8，空間頻率和多普勒頻率超分辨因子ρs和ρd均設(shè)為4，待檢測距離門為300。共選取4種STAP方法進(jìn)行對比，包括傳統(tǒng)的對角加載協(xié)方差求逆法[1,2]，文獻(xiàn)[6]中所提稀疏恢復(fù)方法，文獻(xiàn)[20]中所提聯(lián)合稀疏恢復(fù)算法和本文算法。方便起見，上述4種STAP算法在以下仿真實(shí)驗(yàn)中分別簡稱為 LSMI, SR,SOMP和 SbOMP。需要注意的是，在實(shí)驗(yàn) 1～3中雜噪比(Clutter-to-Noise Ratio, CNR)均設(shè)定為35 dB，實(shí)驗(yàn)1～2中訓(xùn)練樣本數(shù)均選取為12(依照RMB準(zhǔn)則，樣本數(shù)應(yīng)不小于30，因此該情況為樣本嚴(yán)重不足情況)。

4.1 實(shí)驗(yàn)1：雜波空時譜估計(jì)對比

本部分主要對比采用不同算法得到的空時譜分布情況，其中LSMI法采用傳統(tǒng)的MVDR譜。由圖1可以看出，采用稀疏恢復(fù)算法估計(jì)得到雜波譜的分辨率均顯著高于傳統(tǒng) MVDR法。其中，由圖1(b)可以看出，采用 SR法得到的譜分布離散化較為嚴(yán)重，且雜波脊附近區(qū)域存在偽峰；由圖1(c)可以看出，采用SOMP法在雜波脊線上可恢復(fù)得到約r個強(qiáng)雜波點(diǎn)，但在其他區(qū)域也存在大量離散偽峰，這主要是由于噪聲的影響所導(dǎo)致；由圖1(d)可以看出，采用所提SbOMP法可沿雜波脊精確恢復(fù)出連續(xù)的強(qiáng)雜波點(diǎn)，且在其他區(qū)域無偽峰出現(xiàn)。

圖1 采用不同算法估計(jì)空時譜分布情況Fig.1 Space-time spectra corresponding to different estimation algorithm

4.2 實(shí)驗(yàn)2：訓(xùn)練樣本不足情況下雜波抑制性能對比

本節(jié)研究在訓(xùn)練樣本不足情況下上述4種算法的雜波抑制性能對比。本文采用改善因子(Improvement Factor, IF)作為基準(zhǔn)來衡量雜波抑制性能。其中，改善因子定義為輸出信雜噪比(Signal-to-Clutter-plus-Noise Ratio, SCNR)與輸入SCNR的比值，具體表述為：

由圖2可以得出以下3點(diǎn)結(jié)論：一是在訓(xùn)練樣本不足的情況下，采用稀疏恢復(fù)算法的性能要顯著優(yōu)于傳統(tǒng)STAP方法；二是采用聯(lián)合稀疏恢復(fù)算法的主雜波區(qū)的抑制性能要顯著優(yōu)于已有稀疏恢復(fù)算法，這主要是已有稀疏恢復(fù)算法沒有充分利用多個訓(xùn)練樣本中所蘊(yùn)含的雜波信息；三是噪聲環(huán)境下本文所提SbOMB算法的雜波抑制性能要優(yōu)于經(jīng)典聯(lián)合稀疏恢復(fù)算法，這驗(yàn)證了所提算法在噪聲環(huán)境下具有更好的穩(wěn)健性。

圖2 雜波抑制性能對比Fig.2 Clutter suppression performance comparison

4.3 實(shí)驗(yàn)3：不同輸入CNR情況下雜波抑制性能對比

本節(jié)主要研究在不同噪聲背景下，所提基于子空間的聯(lián)合稀疏恢復(fù)算法與經(jīng)典聯(lián)合稀疏恢復(fù)算法的雜波抑制性能對比。其中，Known Cov.算法為已知雜波協(xié)方差矩陣下協(xié)方差矩陣求逆方法，相應(yīng)STAP處理器為最優(yōu)處理器，提供算法性能上限，以方便對比。本實(shí)驗(yàn)中，采用各多普勒通道對應(yīng)改善因子的平均值，即平均改善因子(Average Improvement Factor, AIF)作為標(biāo)準(zhǔn)來衡量雜波抑制性能。共進(jìn)行50次蒙特卡諾仿真。由圖3可以看出，在CNR低于35 dB的情況下，所提SbOMP算法的雜波抑制性能比SOMP算法可改善約3～10 dB；而在噪聲相對較弱(雜噪比較大)情況下，相較于SOMP算法，SbOMP算法的雜波抑制性能也有約1 dB的改善。

5 結(jié)論

本文將經(jīng)典聯(lián)合稀疏恢復(fù)算法應(yīng)用于機(jī)載雷達(dá)STAP處理，并針對其所存在的問題提出了一種基于子空間的聯(lián)合稀疏恢復(fù)STAP方法。研究表明，相比于已有稀疏恢復(fù)算法，采用聯(lián)合稀疏恢復(fù)STAP算法后，對主瓣雜波的抑制性能有了顯著提升；同時，所提基于子空間聯(lián)合稀疏恢復(fù) STAP算法解決了經(jīng)典聯(lián)合稀疏恢復(fù)算法在噪聲背景下性能下降的問題。實(shí)際上，本文算法由于涉及子空間處理，因此運(yùn)算量較經(jīng)典聯(lián)合稀疏恢復(fù)算法有所提升；此外，在極少訓(xùn)練樣本情況下，本文算法的旁瓣抑制性能并不理想，下一步重點(diǎn)針對這些問題開展進(jìn)一步研究。

圖3 不同輸入CNR情況下平均IF對比Fig.3 Average IF comparison of different input CNRs

[1]Klemm R.Principles of Space-Time Adaptive Processing[M].London: Institute of Electical Engineering, 2006.

[2]王永良, 彭應(yīng)寧.空時自適應(yīng)信號處理[M].北京: 清華大學(xué)出版社, 2000.Wang Y L and Peng Y N.Space-Time Adaptive Processing[M].Beijing: Tsinghua University Press, 2000.

[3]Reed I S and Mallett J D.Rapid convergence rate in adaptive arrays[J].IEEE Transactions on Aerospace and Electronics Systems, 1974, 10(6): 853-863.

[4]Selesnick I W, Pillai S U, Li K Y, et al..Angle-Doppler processing using sparse regularization[C]. IEEE International Conference on in Acoustics, Speech and Signal Processing, Dallas, TX, USA, 2010: 2750-2753.

[5]Li J, Zhu X, Stoica P, et al..Iterative space-time adaptive processing[C].IEEE Digital Signal Processing Workshop,Marco Island, FL, 2009: 440-445.

[6]Sun K, Zhang H, Li G, et al..A novel STAP algorithm using sparse recovery technique[C].IEEE International Conference on Geoscience & Remote Sensing Symposium, Cape Town,2009: 336-339.

[7]Sun K, Meng H, Wang Y, et al..Direct data domain STAP using sparse representation of clutter spectrum[J].Signal Processing, 2011, 91(9): 2222-2236.

[8]Sun K, Meng H, Lapierre F D, et al..Registration-based compensation using sparse representation in conformal-array STAP[J].Signal Processing, 2011, 91(10): 2268-2276.

[9]Yang Z C, de Lamare R C, and Li X.l1-Regularized STAP algorithms with a generalized sidelobe canceler architecture for airborne radar[J].IE EE Transactions on Signal Processing, 2012, 60(2): 674-686.

[10]Yang Z C, de Lamare R C, and Li X.Sparsity-aware Space-Time Adaptive Processing algorithms with l1-norm regularization for airborne radar[J].IET Signal Processing,2012, 6(5): 413-423.

[11]Chen S, Donoho D, and Saunders M.Atomic decomposition by basis pursuit[J].SIAM Review, 2001, 43(1): 129-159.

[12]Tropp J A.Just relax: convex programming methods for identifying sparse signals in noise[J].IEEE Transactions on Information Theory, 2006, 53(3): 1030-1051.

[13]Wright S J, Nowak R D, and Figueiredo M A T.Sparse reconstruction by separable approximation[J].IEEE Transactions on Signal Processing, 2009, 57(7): 2479-2493.

[14]Tropp J A.Greed is good: algorithmic results for sparse approximation[J].IEE E Transactions on Information Theory, 2004, 50(10): 2231-2242.

[15]Tropp J A and Gilbert A C.Signal recovery from random measurements via orthogonal matching pursuit[J].IEEE Transactions on Information Theory, 2007, 53(12): 4655-4666.

[16]Needell D and Vershynin R.Signal recovery from incomplete and inaccurate measurements via regularized orthogonal matching pursuit[J].IEEE Journal of Selected Topics on Signal Processing, 2009, 4(2): 310-316.

[17]Gorodnitsky I F and Rao B D.Sparse signal reconstruction from limited data using FOCUSS: a re-weighted minimum norm algorithm[J].IEEE Transactions on Signal Processing,1997, 45(3): 600-616.

[18]Cotter S, Rao B, Engan K, et al..Sparse solutions to linear inverse problems with multiple measurement vectors[J].IEEE Transactions on Signal Processing, 2005, 53(7):2477-2488.

[19]Tropp J, Gilbert A, and Strauss M.Algorithms for simultaneous sparse approximation.Part I: Greedy pursuit[J].Signal Processing, 2006, 86(3): 572-588.

[20]Chen J and Huo X.Theoretical results on sparse representations of multiple-measurement vectors[J].IEEE Transactions on Signal Processing, 2006, 54(12): 4634-4643.

[21]Gribonval R, Rauhut R, Schnass K, et al..Atomsof all channels, unite! average case analysis of multi-channel sparse recovery using greedy algorithms[J].Journal of Fourier Analysis Applications, 2008, 14(5/6): 655-687.

[22]Tropp J A.Algorithms for simultaneous sparse approximation.Part II: Convex relaxation[J].Signal Processing, 2006, 86(3): 589-602.

[23]Wipf D P and Rao B D.An empirical Bayesian strategy for solving the simultaneous sparse approximation problem[J].IEEE Transactions on Signal Processing, 2007, 55(7):3704-3716.

[24]Goldstein J S and Reed I S.Theory of partially adaptive radar[J].IEEE Transactions on Aerospace and Electronic Systems, 1997, 33(4): 1309-1325.