基于幾何方法的結(jié)構(gòu)化協(xié)方差矩陣估計

馬全鑫，杜曉林，董 軍，李建波，田團偉

（1.煙臺大學(xué)計算機與控制工程學(xué)院，山東煙臺 264005；2.重慶郵電大學(xué)通信與信息工程學(xué)院，重慶 400065；3.河南大學(xué)物理與電子學(xué)院，河南開封 475001）

0 引 言

在自適應(yīng)雷達信號處理中，干擾協(xié)方差矩陣（Interference Covariance Matrix，ICM）估計是一個長期存在的基本問題［1-2］。傳統(tǒng)的樣本協(xié)方差矩陣（Sample Covariance Matrix，SCM）估計方法，依賴于均勻環(huán)境下不少于兩倍系統(tǒng)自由度（Degrees of Freedom，DOF）的獨立同分布（Independent and Identically Distributed，IID）訓(xùn)練樣本。然而，在真實場景中雷達所處的環(huán)境通常是非均勻的，以至于僅能獲得較為有限的訓(xùn)練數(shù)據(jù)來估計待檢測單元（Cell Under Test，CUT）的ICM，這會導(dǎo)致估計精度下降，從而影響干擾抑制性能［3］。因此，如何在小樣本情況下準(zhǔn)確估計協(xié)方差矩陣成為當(dāng)前面臨的一個嚴(yán)峻挑戰(zhàn)。

針對上述問題，一類有效的策略是利用雷達場景中的先驗知識提高ICM 的估計精度，其中結(jié)構(gòu)化ICM 估計方法被證實是一種可行的解決方案［4-5］。利用ICM 的結(jié)構(gòu)信息（比如Persymmetric,Toeplitz 結(jié)構(gòu)等），可以減少協(xié)方差矩陣的DOF，從而降低對樣本數(shù)量的依賴性。根據(jù)這一處理范式，國內(nèi)外學(xué)者依據(jù)不同的協(xié)方差模型，提出了多種結(jié)構(gòu)化估計方法［6-13］。文獻［7-9］利用協(xié)方差矩陣的Persymmetric 結(jié)構(gòu)增加了可用訓(xùn)練數(shù)據(jù)，并將其引入廣義似然比（Generalized Likelihood Ratio，GLR）檢測器中，提高了非結(jié)構(gòu)化算法的檢測性能。但上述算法只考慮了矩陣結(jié)構(gòu)信息，當(dāng)結(jié)構(gòu)模型不匹配時可能會導(dǎo)致性能下降。文獻［11］聯(lián)合了Toeplitz 結(jié)構(gòu)信息與雜波環(huán)境知識（比如合成孔徑雷達圖像、數(shù)字高程模型、地理地形圖等），將Toeplitz結(jié)構(gòu)引入知識輔助（Knowledge Aided，KA）色加載矩陣中，提高了算法性能。但該算法依賴于先驗協(xié)方差矩陣的準(zhǔn)確程度，當(dāng)先驗雜波知識失配時會帶來較大的誤差損失。文獻［12］從幾何角度研究了Toeplitz 協(xié)方差矩陣估計問題，證明了當(dāng)訓(xùn)練數(shù)據(jù)較少時，幾何方法能夠獲得優(yōu)越的性能。文獻［14-19］利用額外的約束條件，如正定、低秩和條件數(shù)上限等，以控制所得自適應(yīng)算法的數(shù)值穩(wěn)定性，從而提高了估計精度。

本文遵循幾何范式，提出了基于Persymmetric和Toeplitz 結(jié)構(gòu)的兩種協(xié)方差矩陣估計算法。首先，假設(shè)協(xié)方差矩陣具有上述結(jié)構(gòu)特性，運用這些特性信息對訓(xùn)練數(shù)據(jù)進行處理，進而推導(dǎo)生成了兩種結(jié)構(gòu)樣本協(xié)方差矩陣（Structured Sample Covariance Matrices，SSCMs）。然后，根據(jù)ICM 與相對應(yīng)SSCMs的最小化歐氏距離建立目標(biāo)函數(shù)，并施加正定和條件數(shù)約束。通過該極小化問題的轉(zhuǎn)化，最終求得閉式的協(xié)方差矩陣估計。在分析階段，我們使用兩種場景（空域和多普勒處理）下的輸出信干噪比（Signal-to-Interference-plus-Noise Ratio，SINR）評估了所提出算法的性能。實驗結(jié)果表明，本文所提算法相比于其他同類算法具有更優(yōu)性能。

1 信號模型與結(jié)構(gòu)信息

1.1 信號模型

首先將訓(xùn)練樣本x1,…,xK建模為N維、循環(huán)對稱和零均值的隨機向量，具有相同的協(xié)方差結(jié)構(gòu)

式中，E[·]表示期望，(·)H表示共軛轉(zhuǎn)置，?表示廣義矩陣不等式，M表示ICM，K為樣本總數(shù)。xk∈?N×1的第n個元素表示為xn(k)，則xk定義為

式中，(·)T表示轉(zhuǎn)置。ICM可以具體表示為

式中，Mc表示色干擾矩陣，σ2n表示噪聲功率，IN為N×N維的單位矩陣。

為了使濾波器輸出信號中的SINR 最大化，將最優(yōu)權(quán)矢量定義為［3］

式中，s∈?N×1為目標(biāo)的導(dǎo)向矢量，其表達式取決于具體應(yīng)用場景和雷達配置。在傳統(tǒng)方法中，待估計的協(xié)方差矩陣M由SCM所代替

1.2 結(jié)構(gòu)信息

1.2.1 Persymmetric結(jié)構(gòu)

在雷達系統(tǒng)中，如果采用陣列中心為相位中心，那么ICM 具有Persymmetric 結(jié)構(gòu)，滿足以下等式［6,8］：

式中(·)*表示共軛，JN為N×N維的置換矩陣，即

此外，信號導(dǎo)向矢量s也滿足Persymmetric特性：

1.2.2 Toeplitz結(jié)構(gòu)

與1.2.1 節(jié)相似，線性陣列與均勻發(fā)射脈沖雷達的干擾回波滿足Toeplitz 結(jié)構(gòu)，即協(xié)方差矩陣M沿平行主對角線的每一對角線上的元素都是相同的，此時的ICM可以表示為［6］

可以看出M是由第一行的N個元素所構(gòu)成，我們將其定義為tl,l=0,…,N-1，表示第一行的第l個元素。

2 結(jié)構(gòu)化協(xié)方差矩陣估計

基于Persymmetric 和Toeplitz 結(jié)構(gòu)特性，下文構(gòu)造了兩種SSCMs，并在特定約束集下分別利用它們與ICM 的最小化歐氏距離（也稱為Frobenius 距離）建立優(yōu)化問題，從而提出了兩種估計算法，即Persymmetric 協(xié)方差矩陣估計（Persymmetric Covariance Matrix Estimation,P-CME）算法和Toeplitz 協(xié)方差矩陣估計（Toeplitz Covariance Matrix Estimation,T-CME）算法。

2.1 Persymmetric結(jié)構(gòu)協(xié)方差矩陣估計

基于Persymmetric 結(jié)構(gòu)的協(xié)方差矩陣估計算法可分為兩步：

第一步，根據(jù)協(xié)方差矩陣的Persymmetric 特性，將第k個訓(xùn)練樣本處理為［8］

式中，

是通過分解訓(xùn)練數(shù)據(jù)獲得的獨立樣本向量。此時協(xié)方差矩陣可以表示為

且E[xokxHek]=0。因此，由訓(xùn)練樣本估計的Persymmetric SSCM 為

第二步，由于本文提出的框架依賴于正定矩陣空間中歐式距離的使用，所考慮的矩陣是正定的。因此，必須利用先驗信息得出正定性。為了實現(xiàn)這個目標(biāo)，我們假設(shè)已知噪聲的功率下限σ2（σ2n≥σ2，且為了不失一般性，設(shè)置σ2=0 dB）［14,19］。此外，同時考慮了在條件數(shù)上限約束下最小化M與的Frobenius距離，由此可得優(yōu)化問題

式中，||·||表示Frobenius 范數(shù)，λmax(·)和λmin(·)表示矩陣最大和最小特征值，κM≥1為協(xié)方差矩陣條件數(shù)的上限，可以使用關(guān)于實際雷達電磁環(huán)境的先驗信息或根據(jù)基于觀測的自適應(yīng)框架來指定［12］。

接下來，我們用=UP ΛPUHP表示的譜分解，其中ΛP=diag([d1,d2,…,dN])T，diag(·)表示對角陣，d1≥d2≥… ≥dN是按遞減順序排列的的特征值，UP是一個酉矩陣，其列包含相應(yīng)的特征向量。問題P1的最優(yōu)解為

式中，Λ★=diag([λ★1,λ★2,…,λ★N])T為以下優(yōu)化問題的最優(yōu)解：

式中，Λ=diag([λ1,λ2,…,λN])T，λ1≥λ2≥… ≥λN為M的特征值。為了求解式（16）中的優(yōu)化問題，我們引入輔助變量u＞0，并將問題P＇1進一步等價為

由此可得式（17）的最優(yōu)解為

式中，λ★(u)=[λ1(u),λ2(u),…,λN(u)]T，λi(u)=min(κMu,max(di,max(1,u))),i=1,…,N。此外，μ★為u的最優(yōu)值，可以利用文獻［14］中的結(jié)果以閉式解表示。

2.2 Toeplitz結(jié)構(gòu)協(xié)方差矩陣估計

現(xiàn)在我們考慮第二種情況，即協(xié)方差矩陣M是Toeplitz 結(jié)構(gòu)的。利用Toeplitz 特性和訓(xùn)練樣本將tl估計為［6］

基于估計序列，我們可以構(gòu)造出M的Toeplitz SSCM 估計（用表示）。相似地，通過利用在特定約束集下最小化M與的Frobenius 距離建立優(yōu)化問題，即

由于優(yōu)化問題P1和P2結(jié)構(gòu)的內(nèi)在一致性，故可用P1的求解方法獲得P2的最優(yōu)解。

綜上所述，本文所提出算法的整體流程如下：

步驟1 根據(jù)雷達系統(tǒng)中陣列和脈沖序列先驗知識，判斷ICM 的結(jié)構(gòu)特性（即Persymmetric 或Toeplitz 結(jié)構(gòu)。下述步驟以ICM 滿足Persymmetric結(jié)構(gòu)進行描述）。

步驟2 運用協(xié)方差矩陣的Persymmetric 結(jié)構(gòu)特性并結(jié)合訓(xùn)練數(shù)據(jù)集{x1,…,xK}，生成。

步驟3 通過利用M與的最小化歐氏距離構(gòu)造目標(biāo)函數(shù)||M-||，并遵守MIN和λmax(M)λmin(M) ≤κM約束條件，以此建立優(yōu)化問題P1。

步驟4 求解問題P1，其最優(yōu)解M★=UP Λ★UHP為估計的協(xié)方差矩陣。

步驟5 依據(jù)w★=M★-1s獲得估計的濾波器自適應(yīng)權(quán)矢量。

3 仿真及分析

本節(jié)考慮了兩種典型的雷達信號處理方案：寬帶干擾機干擾接收數(shù)據(jù)的空域處理和雜波干擾回波的多普勒處理［12］。利用SINR 評估了所提出算法的性能，并與現(xiàn)有的一些同類算法進行比較。濾波器的輸出SINR（由200 次蒙特卡洛實驗得到）定義為

式中，|·|表示復(fù)數(shù)的模，=(x)為w的自適應(yīng)估計向量為M的估計量。此外，考慮到實際雷達場景，目標(biāo)狀態(tài)x可以是波達角θ，也可以是歸一化多普勒頻率v。最后，假設(shè)κM=λmax(M)/λmin(M)來進行兩種場景下的模擬實驗［14］。

3.1 空域處理

在此場景下，假定雷達系統(tǒng)配置了一個由N=16 個單元組成的均勻線陣。天線之間的距離為d=λ0/2，其中λ0為波長。假設(shè)總干擾由寬帶干擾機干擾和噪聲組成，則ICM 可以表示為M=Ms+，其中σ2a是噪聲的實際功率水平，Ms是與J個干擾機相關(guān)的協(xié)方差矩陣，即

式中，Bf=B/f0表示相對帶寬，B為所需信號的瞬時帶寬，以及f0=c/λ0，c為光速。此外，δ2i表示第i個干擾機的功率，Φi=2πd(sinθi)λ0表示干擾機相對于天線相位中心的相位角，θi為干擾機的偏離角。那么，在這種情況下產(chǎn)生的導(dǎo)向矢量為s(θ)=[1,ejπsin(θ),…,ejπsin(θ)(N-1)]T。

圖1展示了在高斯訓(xùn)練數(shù)據(jù)下，不同算法SINR 與θ的對應(yīng)關(guān)系。仿真設(shè)置σ2a=0 dB，干擾機數(shù)量J=2，且擁有相同的功率δ2i=20 dB，i=1,2，相位分別為θ1=-30°和θ2=30°，兩個相對帶寬Bf=0.3。圖1（a）與圖1（b）分別將樣本數(shù)目設(shè)為8 和16，并將訓(xùn)練樣本建模為IID、循環(huán)對稱和零均值的高斯隨機向量。圖1將所提算法與文獻中其他算法的性能進行了比較，分別是：基于Frobenius 范數(shù)的估計（Frobenius Norm based Estimator,FNE）和基于譜范數(shù)的估計（Spectral Norm based Estimator,SNE）［14］、約束極大似然估計（Constrained Maximum Likelihood estimator,CML）［16］、快速極大似然估計（Fast Maximum Likelihood estimator,FML）［17］、Oracle 近似收縮估計（Oracle Approximating Shrinkage estimator,OAS）［18］、秩約束極大似然估計（Rank-Constrained Maximum Likelihood estimator,RCML）［19］、Toeplitz 結(jié)構(gòu)約束估計（Toeplitz-Structured Estimator，TSE）［12］以及SCM 算法。由圖1可得，所提算法在SINR 方面都能優(yōu)于同類算法，尤其是在小樣本（K=8）情況下的優(yōu)化效果更為突出。這是因為所提出的算法不僅利用了協(xié)方差矩陣的結(jié)構(gòu)信息，并同時考慮了矩陣的正定特性和條件數(shù)約束。這相當(dāng)于將估計的協(xié)方差矩陣強制約束為趨近于真實協(xié)方差矩陣條件的良好特性，使得未知參數(shù)的不確定區(qū)域減小，以提高估計精度。

圖1 高斯訓(xùn)練數(shù)據(jù)下SINR與θ的對應(yīng)關(guān)系

圖2展示了在非高斯訓(xùn)練數(shù)據(jù)下，不同算法SINR與θ的對應(yīng)關(guān)系。此時的訓(xùn)練樣本建模為

式中，nk～CN(0,σ2aIN)和rk～CN(0,Ms)表示獨立的隨機向量，τk～Γ(0.5,2)為隨機變量。為了進行比較，除了圖1所對比的算法之外，還考慮了兩種用于復(fù)合高斯雜波的算法，即歸一化樣本協(xié)方差矩陣（Normalized Sample Covariance Matrix，NSCM）算法［20］和不動點估計（Fixed-Point Estimator,FPE）算法［21］。由圖2可以看出，所提出的算法在復(fù)合高斯情況下的性能仍優(yōu)于所對比算法。此外，相較于高斯訓(xùn)練數(shù)據(jù)情況，所提算法比同類算法的SINR 增益更高，這表明所提算法在訓(xùn)練數(shù)據(jù)偏離高斯模型的情況下具有較強的穩(wěn)健性。

3.2 多普勒處理

仿真設(shè)置了一個以固定脈沖重復(fù)時間（Pulse Repetition Time,PRT）發(fā)射N=16個相干脈沖序列的雷達系統(tǒng)。在該場景中，x指的是目標(biāo)的歸一化多普勒頻率v=[-1/2,1/2]，導(dǎo)向矢量s(v)=[1,ej2πv,…,ej2πv(N-1)]T。假定雷達是在由海上和地面雜波組成的雙模雜波中工作，則ICM 可以表示為M=Mt+σ2aIN，其中Mt為

式中，fS表示海雜波的歸一化多普勒頻率，CNRS和CNRG分別表示海雜波和地雜波的功率，ρS和ρG分別為海雜波和地雜波的相關(guān)系數(shù)。

圖3展示了在高斯訓(xùn)練數(shù)據(jù)下，不同算法（與圖1中相同）SINR 與v的對應(yīng)關(guān)系。仿真設(shè)置σ2a=10 dB，CNRS=10 dB，CNRG=25 dB，ρS=0.8，ρG=0.95，fS=0.2。訓(xùn)練數(shù)據(jù)滿足高斯分布，且樣本數(shù)目和上一場景設(shè)置相同。觀察曲線可知，在所考慮的情況下，所提出的T-CME 算法在SINR 方面均優(yōu)于對比算法，但P-CME 算法的SINR 略低于TSE算法。然而值得注意的是，Toeplitz 結(jié)構(gòu)為Persymmetric 結(jié)構(gòu)的一個子集［6］，當(dāng)ICM 僅滿足Persymmetric 結(jié)構(gòu)而非Toeplitz 結(jié)構(gòu)時，T-CME 和TSE 算法將不再適用，而本文所提P-CME 算法則優(yōu)于其他同類算法。

圖3 高斯訓(xùn)練數(shù)據(jù)下SINR與v的對應(yīng)關(guān)系

3.3 實測數(shù)據(jù)驗證

本小節(jié)，我們利用實測數(shù)據(jù)進行算法驗證。該實測數(shù)據(jù)來自于2021年海軍航空大學(xué)在煙臺第一海水浴場進行的雷達對海探測數(shù)據(jù)集20210106155330_01_staring.mat，海雜波數(shù)據(jù)具體說明可詳見文獻［22-23］。該雷達系統(tǒng)以固定PRT發(fā)射相干脈沖序列進行對海探測實驗。

我們選取第700 個距離單元作為CUT，使用兩側(cè)臨近距離單元數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)進行實驗。圖4展示了脈沖數(shù)N=16 時，實測數(shù)據(jù)下不同算法SINR 與v的對應(yīng)關(guān)系。由圖4可知，使用實測數(shù)據(jù)的T-CME 算法同樣優(yōu)于同類算法，而P-CME 的SINR 略低于TSE 算法，實測數(shù)據(jù)仿真結(jié)果和模擬數(shù)據(jù)結(jié)果保持一致。

圖4 實測數(shù)據(jù)下SINR與v的對應(yīng)關(guān)系

3.4 復(fù)雜度分析及運行時間對比

計算復(fù)雜度是衡量算法運行效率的重要指標(biāo)，本文所提算法的計算量主要集中在構(gòu)造SSCMs、譜分解及求解u的過程中，通過分析可得出兩種算法生成SSCMs 和進行譜分解的復(fù)雜度為O(N2K+N3+N2)，而u則利用了MATLAB 中一維優(yōu)化函數(shù)fminbnd 進行求解，故所提兩種算法最終的計算復(fù)雜度約為O(N2K+N3+N2)。

表1給出了200 次蒙特卡洛實驗下不同算法的平均運行時間（計算機處理器為Intel Core i5-7200U@2.5 GHz，內(nèi)存為8 GB），由表1可知，在相同仿真條件（采用高斯訓(xùn)練數(shù)據(jù)的空域處理）下，本文所提算法與其他結(jié)構(gòu)化方法運行時間處于同一數(shù)量級，且遠遠小于迭代型方法。

表1 算法運行時間對比結(jié)果（N=16）

3.5 適用性及工程實現(xiàn)分析

本文考慮到在雷達系統(tǒng)中，若線性陣列與脈沖序列滿足對稱或均勻間隔，則產(chǎn)生的ICM 為Persymmetric 或Toeplitz 結(jié)構(gòu)，且合理的利用結(jié)構(gòu)信息可提高協(xié)方差矩陣估計精度。此外，通過上文對算法SINR 和運行時間的對比分析，可得出TCME 算法要優(yōu)于P-CME 算法。因此，在確定ICM具有Toeplitz 結(jié)構(gòu)的性質(zhì)時，優(yōu)先使用T-CME 算法。而當(dāng)ICM 僅滿足Persymmetric 結(jié)構(gòu)時，即可使用P-CME算法。

對于工程實現(xiàn)方面的問題，通過在兩種場景下不同算法性能及運行時間的比較，可以看出所提算法不僅擁有更優(yōu)的干擾抑制能力且計算復(fù)雜度較低，從而有利于信號的實時處理，適合工程應(yīng)用。并且，在實際應(yīng)用中獲取的訓(xùn)練數(shù)據(jù)無論是否服從高斯分布，所提算法都能夠有效地提高雷達干擾抑制性能。

4 結(jié)束語

遵循幾何范式，本文提出了兩種雷達自適應(yīng)信號處理中結(jié)構(gòu)化ICM 估計算法。首先利用結(jié)構(gòu)先驗知識并結(jié)合訓(xùn)練數(shù)據(jù)，構(gòu)建了兩種SSCMs，然后在特定（正定矩陣空間和條件數(shù)上限）約束下通過SSCMs 與ICM 之間的最小化歐氏距離建立優(yōu)化問題并進行求解。所建立的優(yōu)化問題等同于在Frobenius 范數(shù)下將兩種SSCMs 投影到具有實際相關(guān)性的約束集中，以減少估計問題的DOF。在兩種場景下的仿真結(jié)果表明了所提出的算法具有良好的干擾抑制性能。