基于K-均值聚類的SVD雜波抑制算法

黃鳳青，鄭 霖，楊 超，劉爭紅，鄧小芳，扶 明

(桂林電子科技大學廣西無線寬帶通信與信號處理重點實驗室，廣西桂林 541004)

0 引言

強雜波環(huán)境下慢速運動目標的檢測一直是雷達探測的重點之一[1]。由于雷達實際工作環(huán)境比較復雜，目標回波信號易受到強地雜波干擾，使得運動目標很容易被淹沒在雜波中難以檢測[2]。此外，由于運動目標的速度低，目標回波多普勒頻率處在零頻附近，與強地物固定雜波及慢速雜波存在嚴重交疊[3]，傳統(tǒng)的基于動目標顯示(Moving Target Indication,MTI)和動目標檢測(Moving Target Detection，MTD)的雜波抑制算法性能嚴重下降，無法完成對雜波的有效抑制[4]。因此，如何在強雜波背景下實現(xiàn)雜波有效抑制并提高對慢速運動目標的檢測能力成為了一個急需解決的問題[5]。

傳統(tǒng)的雜波抑制方法主要從時域或頻域上濾除雜波，當它們的多普勒頻譜明顯重疊時，不能充分地分離目標回波和雜波。近年來，基于子空間的雜波抑制方法被廣泛應用于雷達目標檢測領域中，現(xiàn)有的子空間方法主要基于奇異值分解(Singular Value Decomposition,SVD)或特征值分解(Eigenvalue Decomposition,EVD)，主要思想是基于回波信號特征，通過特征分解分離目標子空間和雜波子空間，以抑制雜波成分并突顯目標信號。針對海雜波環(huán)境下微小目標檢測問題，文獻[6]提出了正交投影抑制海雜波的方法，通過利用參考單元觀測數(shù)據(jù)構造雜波子空間，然后將待測信號投影到雜波正交子空間中抑制海雜波，能夠降低計算復雜度。針對強雜波環(huán)境下微小目標的檢測問題，文獻[7]以奇異值分解理論為基礎，利用奇異值差分譜的極值點對雷達回波進行分區(qū)，從而達到雜波抑制和目標檢測的目的，文獻[8]提出了一種基于相位編碼及子空間投影的方法，先通過相位解碼白化雜波，再依據(jù)白化雜波和目標的自相關性差異及特征值比值來分離雜波和目標子空間，最后通過子空間投影來抑制雜波。針對血流成像中雜波的抑制問題，文獻[9]對比分析了基于奇異值，空間奇異向量和時間奇異向量選擇閾值的方法，證明了利用空間奇異向量的相關性進行奇異值分解雜波抑制的效果最優(yōu)?；谧涌臻g雜波濾波方法具有一定的雜波抑制效果。然而，現(xiàn)有的基于子空間雜波抑制方法的一個主要問題是雜波子空間難以精確劃分，尤其在慢動目標檢測中，目標回波與雜波存在嚴重混疊，其雜波子空間的準確劃分將會直接影響到目標檢測的性能。如上述在基于子空間方法的濾波中，實現(xiàn)雜波基閾值選擇的方法很多，但是大多數(shù)方法都是根據(jù)特定場景的經驗值來選取雜波基，其并不具有普適性。

針對雜波閾值難以確定的問題，本文設計了一種基于K-均值聚類的子空間投影雜波抑制方法，主要利用奇異值幅度分布、回波信號反映在奇異向量上的空間相關性特征和平均多普勒頻率特征，對相應的特征分量進行聚類。根據(jù)雜波和目標信號具有不同的空時特性，應用 K-均值聚類算法自適應地確認與雜波對應的特征分量，并構建雜波子空間，最后通過正交子空間投影，實現(xiàn)對回波信號中雜波成分的有效抑制，從而提高在強雜波環(huán)境下慢速運動目標的檢測性能。

1 信號模型

雷達發(fā)射機發(fā)射信號為線性調頻信號，雷達脈沖寬度為Tp，脈沖重復周期為Tr，調頻帶寬為B，調頻斜率定義為μ=B/Tp，載頻為fc，則雷達在第n個發(fā)射周期的發(fā)射信號為

(1)

設t=0 時刻，目標相對于雷達的距離為R0，勻速運動的速度為v，則此時目標的回波延遲可以表示為τ=2(R0-vtn)/c,tn=nTr，因此，可以得到單個目標在第n個周期內的回波模型為

exp[j2π(fc(t-nTr-τ)+

(2)

式中，c(t)為雜波，n(t)為高斯白噪聲。

將得到的回波信號xnr(t)重排成快慢時間維回波矩陣，并對快時間維脈沖壓縮可得

(3)

式中，Xc為N×M維的回波數(shù)據(jù)矩陣，N為脈沖積累個數(shù)，M為距離單元數(shù)。

2 基于K-均值聚類的SVD雜波抑制算法

2.1 SVD子空間投影

對經過脈沖壓縮的回波數(shù)據(jù)矩陣Xc進行奇異值分解可得

(4)

式中，C，T，N分別表示雜波、目標和噪聲子空間的奇異值集合，上標H表示矩陣的共軛轉置，ui和vi分別為矩陣Xc的第i個左奇異值向量和右奇異值向量，λi為第i個奇異值，按從大到小降序排列。在確定表征雜波的奇異值集合C之后，可以得到下面的奇異向量矢量積來表示雜波子空間：

(5)

(6)

式中，I為單位矩陣。為了抑制強雜波的影響，將回波信號矩陣Xc投影到該正交子空間中：

(7)

2.2 雜波子空間的確定

為了充分抑制雜波，如何選取表征雜波子空間所對應的奇異譜向量是關鍵問題。現(xiàn)有奇異值分解大多數(shù)只從奇異值分布情況(如最大奇異值或奇異值差分譜)來確定雜波基，或單從奇異向量入手確定雜波基，未能充分利用雜波的特性來選取雜波基，導致雜波抑制出現(xiàn)較大的殘余。本文分析了回波數(shù)據(jù)中主要構成成分的特性，接收端回波信號主要由雜波、目標信號以及噪聲三種成分構成，其中雜波成分是由發(fā)射信號經靜止物體如地面、墻體、灌木等或慢速運動的物體一次或多次反射造成，主要分布在零多普勒附近，因此在不同的距離單元具有較強的空間相關性且平均多普勒較低，同時由于墻體和地面的反射能力強，雜波的回波能量強度較高；而待檢測運動目標只在某個或若干個距離單元出現(xiàn)且目標的RCS比較小，所以目標的空間相關性較小、回波能量弱，且平均多普勒一般高于靜止或慢速雜波信號。根據(jù)雜波和信號的差異性，利用奇異值幅值分布、奇異向量的空間相關性和平均多普勒頻率三個特性對特征分量進行聚類，無需人為設定閾值參數(shù)估計雜波基，能夠自適應確定雜波子空間所對應的奇異向量。

2.2.1 奇異值譜及奇異向量分析

1) 奇異值譜分布。在強雜波環(huán)境下，雜波占主導地位，雜波回波能量往往比目標信號回波能量大，通過奇異值分解獲得雷達回波數(shù)據(jù)矩陣的歸一化奇異值幅度分布如圖1(a)所示，雖然只截取了前40個奇異值，但可以看到奇異值的幅值下降得比較平緩，不存在明顯的突變點，這是因為雜波成分的特征分布具有多維性，導致雜波子空間和目標子空間發(fā)生混疊。

2) 奇異向量的相關性。回波矩陣Xc的奇異值分解結果是V奇異向量提供空間信息，而U奇異向量表示Xc中的多普勒信息[10-11]。地雜波的相關性主要是指距離項的空間相關性[12]，不同距離分辨單元中的地雜波具有較強的相關性，而目標只存在某一個或若干個距離分辨單元中，因此目標信號的空間相關性較差，利用目標信號分量和地雜波分量的空間相關性差異，也是子空間分解類雜波抑制算法能夠有效抑制雜波的基礎[13]。接收的回波信號中雜波與雜波之間的相關性較強，而雜波與目標信號回波的相關性較弱，對回波信號矩陣進行奇異值分解，得到的奇異向量表征了這種特性。其中，因為回波信號中雜波成分占主導地位，雜波能量較強，因此最大奇異值對應的奇異向量V1表征了雜波成分，計算V1與左奇異向量之間的相關系數(shù)：

(8)

式中：Vk為第k個奇異值對應的右奇異向量，k∈[1,M]，M為快時間維的采樣點數(shù)；V1為第一個奇異值對應的右奇異向量；E(Vk)，σk分別表示Vk的均值和標準差;ρV1,Vk為相關系數(shù)向量。相關系數(shù)向量的曲線圖如圖1(b)所示，從圖中可知，前面二十幾個奇異向量的相關性大于后面奇異向量的相關性，揭示了與雜波所對應的低階奇異向量的高相關性，以及對應目標的高階奇異向量的低空間相關性。

3) 平均多普勒估計。地雜波主要分布在零多普勒頻率附近，相對于運動目標，平均多普勒較小。利用滯后一位的自相關估計器，計算每個左奇異向量Uk的平均多普勒頻率：

(9)

(10)

圖1 奇異值及奇異向量特征分布圖

2.2.2K-均值聚類

K-均值聚類算法主要通過迭代法將一個給定的數(shù)據(jù)集劃分為k個簇，具有簡單高效、收斂速度快等優(yōu)點，被廣泛應用于科學研究等諸多領域。K-均值算法基本思想是首先確定期望被劃分簇的個數(shù)k值，然后從給定的數(shù)據(jù)集中隨機選取k個數(shù)據(jù)點作為質心即初始聚類中心，一般以歐式距離作為相似性度量計算數(shù)據(jù)集每一個數(shù)據(jù)點到k個質心的距離，最后根據(jù)最近鄰原則將數(shù)據(jù)點劃歸到就近質心的簇中。接下來重新確定k個簇中的質心，計算式如式(11)，假設有數(shù)據(jù)集X，將其按照上述過程劃分為k個簇(C1,C2,…,Ck),則求質心μi的表達式為

(11)

將新得到的質心與上一次劃分簇時的質心相比較，判定質心是否發(fā)生變化，如果沒有發(fā)生變化，則迭代終止得到最終的聚類質心，否則根據(jù)上述方法重新劃分簇，直到滿足算法的收斂條件為止。

依據(jù)奇異值幅度分布、奇異向量的空間相關性特征和平均多普勒頻率特征三個特征，本文采用K-均值聚類算法對回波信號奇異值分解后的數(shù)據(jù)進行聚類，如圖2所示。從圖中可知，“聚類3”數(shù)據(jù)(藍色)分布在奇異值幅值較大，空間相關性較強，以及平均多普勒較小的區(qū)域，與雜波的特性相符，因此“聚類3”被確定為代表雜波成分的類。而“聚類1”和“聚類2”的數(shù)據(jù)分布在奇異值幅值較小、空間相關性較弱和平均多普勒較大的區(qū)域，因此分別代表了目標回波成分與噪聲成分。因為噪聲能量較小，且噪聲的頻率豐富，在整個頻譜上都有成分，而弱目標的回波信號能量小，目標的多普勒信息未知，因此聚類為表征目標回波與噪聲的類的類間距非常接近，不易分開。所以本文選擇表征為雜波信號的聚類作為雜波子空間,如圖3所示，藍色部分奇異值對應的奇異向量作為雜波基。基于此，本文自適應地確定代表雜波成分的奇異向量作為雜波基，然后利用構建的雜波基，根據(jù)最小二乘法做子空間投影以實現(xiàn)抑制回波信號中的雜波成分。

圖2 K-均值聚類三維分布圖

圖3 奇異值分布圖

3 仿真實驗

為了驗證算法的有效性，本文按照表1所給的雷達參數(shù)進行了仿真實驗，并與現(xiàn)有的基于傳統(tǒng)雜波抑制方法進行了比較分析。仿真實驗的雜波信號采用符合地雜波特性的瑞利雜波模型生成，且預先在場景中設置一個徑向速度為2 m/s,距離為40 m運動目標。為了體現(xiàn)強雜波環(huán)境，設置的信雜比為-28 dB。

表1 雷達系統(tǒng)仿真參數(shù)

圖4(a)給出了僅相參積累后的距離-速度圖。從圖中可以看出，由于地雜波的擴展使得目標信號在多普勒譜上被雜波成分所掩蓋，目標無法被清晰地檢測出來。圖4(b)、(c)、(d)分別是采用二階MTI算法、傳統(tǒng)子空間方法和本文基于K-均值聚類的子空間雜波抑制方法處理結果圖。MTI處理的基本思想是基于雜波與目標回波多普勒頻率的不同，在零頻附近產生一個凹陷來濾除零頻雜波，但對于慢速目標來說，其回波多普勒也處在零頻附近，MTI算法在抑制零頻雜波的同時，對目標回波也有較大的衰減。而且在強雜波環(huán)境下目標被雜波所掩沒時，MTI處理后的殘余雜波對運動目標仍有較大影響，如圖4(b)所示，MTI算法無法起到較好的抑制效果。傳統(tǒng)的子空間方法一般選擇最大的一個或幾個特征值對應的特征向量來構建雜波子空間，但是雜波子空間往往存在交織特性，當奇異值閾值選擇過小時，抑制后仍有較大雜波殘留；當奇異值閾值選擇過大時，部分目標信號的成分會被當作雜波成分抑制掉，導致目標能量減弱，不利于目標檢測。如圖4(c)所示選取前15個特征向量作為雜波基，經過子空間投影后仍殘留較多的雜波，雜波抑制效果并不理想。

圖4(d)是本文的子空間方法雜波抑制的效果，從圖上可以看出目標的距離為40 m,速度為2 m/s，和預先設置的目標運動參數(shù)一致，表明雜波能夠得到較好的雜波抑制效果且目標信號得到凸顯。本文的子空間方法主要是從回波信號特征出發(fā)，利用奇異值分解后的奇異值和奇異向量特性，結合聚類算法自適應地選取與雜波所對應的特征成分，然后通過構建雜波子空間并進行子空間投影，實現(xiàn)對雜波的有效抑制。

(a) 未進行雜波抑制處理

(b) MTI

(c) 傳統(tǒng)子空間方法

(d) 本文SVD法圖4 多種目標檢測算法抗雜波的性能比較

為了驗證算法在多目標場景下的適用性，在實驗中預先設置徑向速度為-2 m/s、距離為30 m和速度為2 m/s、距離為30 m的兩個運動目標。如圖5(a)所示，在本文場景下運動目標均被雜波所掩蓋。

在圖5(b)中，MTI方法在抑制零頻雜波的同時，對運動目標回波也有較大衰減。在圖5(c)中經過傳統(tǒng)子空間方法處理后兩個目標得以凸顯，但是由于雜波基選取不恰當，仍存在較多的雜波殘留。而圖5(d)本文子空間方法可以看出，基于K-均值聚類來自適應選取奇異向量作為雜波基，能夠更加充分地抑制掉雜波，使目標能更加清晰地被檢測出來。

(a) 未經過雜波抑制處理

(b) MTI方法

(c) 傳統(tǒng)子空間方法

(d) 本文子空間方法圖5 多目標場景下目標檢測算法抗雜波的性能比較

4 結束語

針對強雜波環(huán)境下慢速運動目標被掩蓋的問題，本文設計了一種基于K-均值聚類的子空間投影雜波抑制算法。相比于傳統(tǒng)子空間方法只選取最大一個或幾個特征值所對應的特征向量構建雜波子空間，本文算法依據(jù)回波信號特性，得到相應的奇異值譜分布，以及奇異向量的空間相關性特征和平均多普勒頻率特征三個特征，采用K-均值聚類算法對奇異值分解后的數(shù)據(jù)進行聚類，自適應地確定相應的奇異向量來構建雜波子空間。本文方法無需人為設定閾值參數(shù)來估計雜波基，實現(xiàn)了雜波的自適應抑制，也為在雜波基的選取上提供了一種新的思路。實驗結果也表明了相對于傳統(tǒng)的雜波抑制方法，本文方法提高了強雜波環(huán)境下對慢速運動目標的檢測性能。