基于幾何方法的結構化協方差矩陣估計

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

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

0 引 言

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

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

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

1 信號模型與結構信息

1.1 信號模型

首先將訓練樣本x1,…,xK建模為N維、循環對稱和零均值的隨機向量，具有相同的協方差結構

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

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

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

為了使濾波器輸出信號中的SINR 最大化，將最優權矢量定義為［3］

式中，s∈?N×1為目標的導向矢量，其表達式取決于具體應用場景和雷達配置。在傳統方法中，待估計的協方差矩陣M由SCM所代替

1.2 結構信息

1.2.1 Persymmetric結構

在雷達系統中，如果采用陣列中心為相位中心，那么ICM 具有Persymmetric 結構，滿足以下等式［6,8］：

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

此外，信號導向矢量s也滿足Persymmetric特性：

1.2.2 Toeplitz結構

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

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

2 結構化協方差矩陣估計

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

2.1 Persymmetric結構協方差矩陣估計

基于Persymmetric 結構的協方差矩陣估計算法可分為兩步：

第一步，根據協方差矩陣的Persymmetric 特性，將第k個訓練樣本處理為［8］

式中，

是通過分解訓練數據獲得的獨立樣本向量。此時協方差矩陣可以表示為

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

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

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

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

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

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

由此可得式（17）的最優解為

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

2.2 Toeplitz結構協方差矩陣估計

現在我們考慮第二種情況，即協方差矩陣M是Toeplitz 結構的。利用Toeplitz 特性和訓練樣本將tl估計為［6］

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

由于優化問題P1和P2結構的內在一致性，故可用P1的求解方法獲得P2的最優解。

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

步驟1 根據雷達系統中陣列和脈沖序列先驗知識，判斷ICM 的結構特性（即Persymmetric 或Toeplitz 結構。下述步驟以ICM 滿足Persymmetric結構進行描述）。

步驟2 運用協方差矩陣的Persymmetric 結構特性并結合訓練數據集{x1,…,xK}，生成。

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

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

步驟5 依據w★=M★-1s獲得估計的濾波器自適應權矢量。

3 仿真及分析

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

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

3.1 空域處理

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

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

圖1展示了在高斯訓練數據下，不同算法SINR 與θ的對應關系。仿真設置σ2a=0 dB，干擾機數量J=2，且擁有相同的功率δ2i=20 dB，i=1,2，相位分別為θ1=-30°和θ2=30°，兩個相對帶寬Bf=0.3。圖1（a）與圖1（b）分別將樣本數目設為8 和16，并將訓練樣本建模為IID、循環對稱和零均值的高斯隨機向量。圖1將所提算法與文獻中其他算法的性能進行了比較，分別是：基于Frobenius 范數的估計（Frobenius Norm based Estimator,FNE）和基于譜范數的估計（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 結構約束估計（Toeplitz-Structured Estimator，TSE）［12］以及SCM 算法。由圖1可得，所提算法在SINR 方面都能優于同類算法，尤其是在小樣本（K=8）情況下的優化效果更為突出。這是因為所提出的算法不僅利用了協方差矩陣的結構信息，并同時考慮了矩陣的正定特性和條件數約束。這相當于將估計的協方差矩陣強制約束為趨近于真實協方差矩陣條件的良好特性，使得未知參數的不確定區域減小，以提高估計精度。

圖1 高斯訓練數據下SINR與θ的對應關系

圖2展示了在非高斯訓練數據下，不同算法SINR與θ的對應關系。此時的訓練樣本建模為

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

3.2 多普勒處理

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

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

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

圖3 高斯訓練數據下SINR與v的對應關系

3.3 實測數據驗證

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

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

圖4 實測數據下SINR與v的對應關系

3.4 復雜度分析及運行時間對比

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

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

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

3.5 適用性及工程實現分析

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

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

4 結束語

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