溫 柔，孫惜媛，何 姿，顧鵬飛，丁大志

(南京理工大學(xué) 電子工程與光電技術(shù)學(xué)院，南京 210094)

0 引言

雷達(dá)信號(hào)的目標(biāo)識(shí)別技術(shù)在軍事檢測(cè)防御系統(tǒng)中意義非凡，盲源分離(blind source separation, BSS)是解決雷達(dá)系統(tǒng)中混疊信號(hào)的有效途徑[1]。BSS最早是從“雞尾酒會(huì)問(wèn)題”引申而來(lái)，其本質(zhì)上是一種功能強(qiáng)大的信號(hào)處理問(wèn)題。該問(wèn)題描述為：在若干個(gè)麥克風(fēng)同時(shí)接收到若干個(gè)聲源的情況下，從混合的麥克風(fēng)信號(hào)中恢復(fù)出各個(gè)聲源，欠定盲源分離(underdetermined blind source separation, UBSS)可以理解為麥克風(fēng)數(shù)量小于聲源數(shù)量。本文在雷達(dá)系統(tǒng)的背景下，從接收到的微動(dòng)群目標(biāo)的混合信號(hào)中分離出各目標(biāo)的回波信號(hào)。

自1985年Herault和Jutten[2]首次提出至今，BSS一直是研究的熱門(mén)方向，該領(lǐng)域中應(yīng)用最廣泛的是Comon提出的獨(dú)立成分分析(independent component analysis, ICA)方法[3]。由于UBSS的不可逆性， ICA方法并不適用于此場(chǎng)景，因此Lewicki根據(jù)源信號(hào)的稀疏特性提出了稀疏分量分析方法(sparse component analysis, SCA)[4]，成為欠定領(lǐng)域中的前沿方法并被廣泛應(yīng)用[5-6]。

2001年，Bofill和Zibulevsky基于SCA提出了經(jīng)典的“兩步法”[7]，該方法的許多參數(shù)設(shè)置雖然在當(dāng)時(shí)沒(méi)有具體的理論支持，但仍具有重要的貢獻(xiàn)和參考價(jià)值。一般來(lái)說(shuō)，兩步法首先通過(guò)未知的混合矩陣將若干個(gè)源信號(hào)混合得到回波信號(hào)，在該回波信號(hào)的基礎(chǔ)上對(duì)混合矩陣進(jìn)行估計(jì)，然后利用前面估計(jì)出來(lái)的混合矩陣分離出各源信號(hào)。在第一步中，文獻(xiàn)[8]利用混合信號(hào)的局部協(xié)方差矩序列的非平穩(wěn)特性，將混合矩陣估計(jì)轉(zhuǎn)化為聯(lián)合近似對(duì)角化(joint approximate diagonalization, JAD) 最小化模型，文獻(xiàn)[9]將JAD最小化問(wèn)題轉(zhuǎn)化為張量分解問(wèn)題。文獻(xiàn)[10-11]使用合適的算法來(lái)檢測(cè)單源點(diǎn)，利用稀疏的單源點(diǎn)并結(jié)合聚類(lèi)算法進(jìn)行混合矩陣估計(jì)。文獻(xiàn)[12]提出了經(jīng)典的K-means聚類(lèi)方法用于估計(jì)混合矩陣，該方法主要的不足之處在于源信號(hào)的個(gè)數(shù)要已知，這并不符合實(shí)際情況。由于欠定情況下的解混操作沒(méi)有唯一解，第二步將源分離問(wèn)題建立為各種優(yōu)化模型。最經(jīng)典的方法為線性規(guī)劃算法中的l1范數(shù)最小化方法和最短路徑法[13-14]，l1范數(shù)最小化方法要求源信號(hào)足夠稀疏，最短路徑法只適用于在兩個(gè)觀測(cè)通道的情況下，因此這些經(jīng)典方法都具有非常大的局限性。文獻(xiàn)[15]提出了一種改進(jìn)的l1范數(shù)最小化方法，以最小化源信號(hào)中分量值較小的部分元素為目標(biāo)函數(shù)，采用線性規(guī)劃方法求解新的數(shù)學(xué)優(yōu)化模型的精確解，實(shí)現(xiàn)源信號(hào)重構(gòu)。文獻(xiàn)[16]基于源的拉普拉斯分布假設(shè)，提出了一種基于lp范數(shù)的源分離優(yōu)化模型，但是，當(dāng)p小于0.75時(shí)，求得的解為局部極小值而非全局最優(yōu)解。第二步的另一種經(jīng)典方法是壓縮感知理論，文獻(xiàn)[17]在用壓縮感知理論解決欠定盲源分離問(wèn)題時(shí)提出了一種基于分層耦合的字典訓(xùn)練方法，通過(guò)對(duì)先驗(yàn)信息的分層訓(xùn)練，實(shí)現(xiàn)雷達(dá)信號(hào)的高精度重構(gòu)。

為解決上述問(wèn)題，本文在“兩步法”框架的基礎(chǔ)上，以空中微動(dòng)群目標(biāo)為源信號(hào)[18]，雷達(dá)接收的回波為混合信號(hào)，提出了一種改進(jìn)的兩階段混合矩陣估計(jì)和源分離方法。第一步，利用混合信號(hào)的局部協(xié)方差序列，將混合矩陣估計(jì)問(wèn)題轉(zhuǎn)化為特征向量聚類(lèi)問(wèn)題。在聚類(lèi)中，通過(guò)檢查數(shù)據(jù)集中每點(diǎn)的R鄰域來(lái)實(shí)現(xiàn)簇的分類(lèi)。如果點(diǎn)I的R鄰域中數(shù)據(jù)點(diǎn)的個(gè)數(shù)多于P，即I點(diǎn)周?chē)芏容^大時(shí)，則創(chuàng)建一個(gè)以I為聚類(lèi)中心的簇，否則I點(diǎn)視為離群點(diǎn)，最后通過(guò)求得每一簇的均值來(lái)實(shí)現(xiàn)混合矩陣的估計(jì)。該方法優(yōu)點(diǎn)在于不需要預(yù)先輸入簇的個(gè)數(shù)，并能有效處理離群點(diǎn)。第二步，將源分離問(wèn)題轉(zhuǎn)化為基于lp范數(shù)(0

1 盲源分離系統(tǒng)模型

BSS的系統(tǒng)模型可以劃分為線性和非線性?xún)煞N混合模型[19]，每種模型又分為瞬時(shí)、延時(shí)和卷積混合模型。本文主要研究的是線性瞬時(shí)混合模型。

在不考慮噪聲的情況下，線性瞬時(shí)混合模型的BSS過(guò)程可以表示為：

X(t)=HS(t)

(1)

Y(t)=WX(t)

(2)

式(1)中，S(t)=[s1(t),s2(t),…,sN(t)]T為N維的未知源信號(hào)，H為M*N維的混合矩陣，X(t)=[x1(t),x2(t),…,xM(t)]T為M維的混合信號(hào)。式(2)中，W為N*M維的分離矩陣，Y(t)=[y1(t),y2(t),…,yN(t)]T為N維的分離信號(hào)。

以分離彈道中段的彈頭群目標(biāo)為研究對(duì)象，當(dāng)?shù)孛胬走_(dá)檢測(cè)到多個(gè)目標(biāo)的混疊回波時(shí)，為了獲取目標(biāo)的有用信息，需要從混疊回波中將每一個(gè)目標(biāo)分離出來(lái)。如圖1所示，顯示了空中微動(dòng)群目標(biāo)的BSS模型圖。該模型由M個(gè)雷達(dá)系統(tǒng)接收N個(gè)目標(biāo)的回波，從而獲得M個(gè)通道的混合信號(hào)，然后利用BSS算法分離出每一個(gè)目標(biāo)的回波。

圖1 BSS模型圖Fig.1 The model diagram of BSS

為了解決上述模型的問(wèn)題，本文的方法需要建立在兩個(gè)假設(shè)前提下：

1) 混合信號(hào)的數(shù)量小于源信號(hào)的數(shù)量，也就是說(shuō)需要在欠定的情況下，即M

2) 每個(gè)源信號(hào)需要存在局部?jī)?yōu)勢(shì)[20-21]，即在時(shí)域上至少存在一個(gè)子塊是僅由一個(gè)源信號(hào)起主導(dǎo)作用的。在不同的微動(dòng)參數(shù)，例如微動(dòng)頻率、微動(dòng)形式等，以及不同的尺寸、樣式下，彈頭群目標(biāo)的回波會(huì)存在局部?jī)?yōu)勢(shì)，對(duì)于任意一個(gè)目標(biāo)源信號(hào)僅需保證在某一時(shí)間段內(nèi)占主導(dǎo)作用，即可將該信號(hào)從混合信號(hào)中分離出來(lái)。如圖2所示，在三個(gè)源信號(hào)的電場(chǎng)實(shí)部圖中，紅色框中的每一塊在其時(shí)域中僅存在一個(gè)源信號(hào)有電場(chǎng)分量。

圖2 局部?jī)?yōu)勢(shì)Fig.2 Local dominance

2 算法原理

2.1 整體框圖

本文將用一種改進(jìn)的“兩步法”來(lái)解決圖2所示的BSS模型，該方法的原理框圖如圖3所示。

圖3 兩步法框圖Fig.3 The framework of ‘two-stage’

為了簡(jiǎn)單起見(jiàn)，在下面的討論中X表示混合信號(hào)、S表示源信號(hào)、Y表示分離信號(hào)。在第一步中，先將混合信號(hào)X分塊，構(gòu)造每一塊的局部協(xié)方差矩陣，然后對(duì)每一個(gè)局部協(xié)方差矩陣進(jìn)行特征值分解(eigenvalue decomposition, EVD)，提取出最大特征值對(duì)應(yīng)的特征向量，最后將每一塊提取的特征向量進(jìn)行聚類(lèi)，得到混合矩陣。在第二步中，先對(duì)分離信號(hào)進(jìn)行初始化，然后將混合信號(hào)與估計(jì)的混合矩陣進(jìn)行列式分解，通過(guò)循環(huán)迭代得到分離信號(hào)。本文用密度聚類(lèi)[22]的方法來(lái)完成第一步混合矩陣的估計(jì)，再用迭代拉格朗日乘子法[23]來(lái)完成第二步源信號(hào)的分離，以下章節(jié)將給出每一步的詳細(xì)講解。

2.2 混合矩陣估計(jì)

2.2.1 數(shù)據(jù)預(yù)處理

首先利用源信號(hào)的非平穩(wěn)特性引入局部協(xié)方差矩陣。將M*L維的混合信號(hào)X平均劃分為Q個(gè)相互獨(dú)立的子塊，每個(gè)子塊相當(dāng)于一個(gè)M*[L/Q]維的矩陣。定義混合信號(hào)第q個(gè)子塊的局部協(xié)方差矩陣為：

(3)

(4)

(5)

對(duì)每一塊的局部協(xié)方差矩陣進(jìn)行EVD為：

(6)

式(6)中，Vq是特征向量矩陣，∑q為特征值矩陣。提取Vq中最大特征值所對(duì)應(yīng)的特征向量為aq，對(duì)所有Q個(gè)子塊進(jìn)行如上所示的EVD，得到Q個(gè)特征向量，組成特征向量矩陣為A=[a1,…,aQ]。

對(duì)矩陣A的各列向量進(jìn)行密度聚類(lèi)，找出N個(gè)簇，將每一簇取均值求得混合矩陣。在(M,N)=(2,3)的情況下，將所有求得的特征向量降維映射到二維空間中，其散點(diǎn)圖如圖4所示。圖中共有3個(gè)簇，不同簇被標(biāo)記成不同的顏色，每一個(gè)簇的均值對(duì)應(yīng)估計(jì)的混合矩陣每個(gè)列向量，不在大密度區(qū)域中的點(diǎn)被標(biāo)記為離群點(diǎn)。

圖4 聚類(lèi)示意圖Fig.4 The diagram of the clustering

2.2.2 密度聚類(lèi)算法

密度聚類(lèi)算法核心思想是用一個(gè)點(diǎn)R鄰域內(nèi)的鄰點(diǎn)數(shù)量來(lái)衡量該點(diǎn)所在空間的密度，將大密度點(diǎn)及周?chē)芏瓤蛇_(dá)的鄰點(diǎn)劃分為一簇，從而找出形狀不規(guī)則的簇，且不需要提前輸入簇的個(gè)數(shù)。

下面給出一些基本定義，對(duì)于特征向量矩陣A的每一列特征向量a，稱(chēng)之為一個(gè)數(shù)據(jù)點(diǎn)。本算法重點(diǎn)考慮了兩個(gè)因素，定義密度時(shí)的鄰域半徑R和定義聚類(lèi)中心點(diǎn)時(shí)的最小鄰域點(diǎn)數(shù)P。

1)R鄰域：定義NR(a)={b|d(b,a)≤R}為a的R鄰域，d(b,a)表示b與a之間的距離，顯然a,b∈A。

2)密度ρ：定義ρ(a)=|NR(a)|為a點(diǎn)的密度，它是依賴(lài)于鄰域半徑R的整數(shù)值。

3)聚類(lèi)中心：若ρ(a)≥P，則稱(chēng)a為聚類(lèi)中心，由聚類(lèi)中心構(gòu)成的集合為AC,ANC為非聚類(lèi)中心點(diǎn)構(gòu)成的集合。

4)離群點(diǎn)：若ρ(a)

5)密度直達(dá)：若a∈AC,b∈NR(a)，則稱(chēng)b是從a出發(fā)密度直達(dá)的。

6)密度可達(dá)：設(shè)a1,a2,…,am∈A，m≥2，若ai+1是從ai密度直達(dá)的，i=1,2,…,m-1，則am是從a1密度可達(dá)的。

密度聚類(lèi)算法的過(guò)程基本描述為：任意滿(mǎn)足密度條件的聚類(lèi)中心點(diǎn)I，數(shù)據(jù)集中所有從I點(diǎn)密度可達(dá)的點(diǎn)O所組成的集合構(gòu)成了一個(gè)完整的簇,接著繼續(xù)選擇另一個(gè)聚類(lèi)中心點(diǎn)去尋找密度可達(dá)的數(shù)據(jù)集，構(gòu)成另一個(gè)簇，直到?jīng)]有新的點(diǎn)加到任何簇中。根據(jù)密度聚類(lèi)算法估計(jì)混合矩陣的具體過(guò)程如下：

步驟一：輸入2.2.1節(jié)中的特征向量矩陣A，設(shè)置鄰域參數(shù)R和P，初始化聚類(lèi)中心集，簇的個(gè)數(shù)。對(duì)于數(shù)據(jù)中的每個(gè)點(diǎn)，計(jì)算其對(duì)應(yīng)的k個(gè)最近鄰域距離(一般k設(shè)置為4)，并將數(shù)據(jù)集中所有點(diǎn)對(duì)應(yīng)的k個(gè)最近鄰域距離按照降序方式排序，繪制成“K-距離”曲線。R的值一般設(shè)定為曲線圖拐點(diǎn)位置的值。P值的設(shè)定需要滿(mǎn)足：P≥dim+1，dim表示數(shù)據(jù)的維度。例如，在本文算例中R設(shè)定為0.03，P設(shè)定為3。

步驟二：對(duì)于A中任意的a，通過(guò)歐式距離的度量方式，找到其R鄰域子集NR(a)，若子集個(gè)數(shù)|NR(a)|≥P，則將a加入聚類(lèi)中心集AC，否則標(biāo)記a為離群點(diǎn)。

步驟三：若AC=φ，則算法結(jié)束，否則轉(zhuǎn)入步驟四。

步驟四：在AC中選擇一個(gè)聚類(lèi)中心a，建立新簇K，將a的R鄰域子集NR(a)所有的點(diǎn)加入K。并對(duì)NR(a)中所有未被處理的點(diǎn)b檢查其R鄰域子集NR(b)，若|NR(b)|≥P，則將NR(b)中未歸入任意簇的對(duì)象加入簇K。

步驟五：輸出N個(gè)簇：K1，…,KN，對(duì)每一個(gè)簇求均值，得到C1,…,CN，則C=[C1,…,CN]為估計(jì)的混合矩陣。

2.3 源信號(hào)分離

本文采用一種基于lp范數(shù)最小化的模型來(lái)實(shí)現(xiàn)源信號(hào)的分離，其目標(biāo)在于根據(jù)線性瞬時(shí)混合模型式(1)找到Sd的最稀疏項(xiàng)，因此源分離的等效模型可以表示為[23]：

(7)

用拉格朗日乘子將式(7)轉(zhuǎn)化為一個(gè)新的無(wú)約束優(yōu)化問(wèn)題：

(8)

式(8)中，α∈M為拉格朗日乘子。式(8)進(jìn)一步展開(kāi)為：

Sd=Ψ-1(Sd)CH(CΨ-1(Sd)CH)-1Xd

(9)

(10)

(11)

3 仿真分析

3.1 評(píng)價(jià)標(biāo)準(zhǔn)

本文通過(guò)相似度和相對(duì)均方根誤差(RRMSE)來(lái)衡量BSS系統(tǒng)的性能。

相似度是用來(lái)衡量?jī)蓚€(gè)信號(hào)波形之間的相關(guān)程度的，表達(dá)式為：

(12)

相似度的范圍是0≤ρi≤1，ρi越大說(shuō)明分離的效果越好，反之則越差，當(dāng)ρi=1時(shí)表示分離的信號(hào)與源信號(hào)完全一樣，一般來(lái)說(shuō)，相似度大于0.8時(shí)則可以認(rèn)為分離效果較好。

相對(duì)均方根誤差(RRMSE)用來(lái)描述兩個(gè)信號(hào)之間的偏離程度，表達(dá)式為：

(13)

RRMSE越小說(shuō)明分離的效果越好，分離出的信號(hào)數(shù)據(jù)與原始信號(hào)數(shù)據(jù)之間相差越小，相反的，分離效果越差。

3.2 仿真算例

本文算例采用球頭錐的模型來(lái)模擬真實(shí)彈頭結(jié)構(gòu)，利用旋轉(zhuǎn)對(duì)稱(chēng)體球面等效源區(qū)域分解算法(EPA-BOR)[24]得到彈頭的回波信號(hào)，分別用本文提出的算法和文獻(xiàn)[25-26]提出的算法進(jìn)行混合矩陣估計(jì)和源信號(hào)的分離。

考慮3個(gè)目標(biāo)兩個(gè)通道的情況，設(shè)目標(biāo)回波中有3個(gè)彈頭目標(biāo)作進(jìn)動(dòng)，每個(gè)彈頭的大小各不相同，彈頭質(zhì)心為目標(biāo)的中心，彈頭目標(biāo)模型如圖5所示。

圖5 彈頭目標(biāo)模型Fig.5 The model of the warhead targets

彈頭的尺寸和微動(dòng)參數(shù)如表1所列。

表1 目標(biāo)幾何和微動(dòng)參數(shù)Tab.1 Target geometry and fretting parameters

雷達(dá)發(fā)射采用單頻脈沖信號(hào)，載頻為10GHz，采樣頻率為500Hz，雷達(dá)視線角為70°，觀測(cè)時(shí)間為 4s。一般來(lái)說(shuō)，目標(biāo)相對(duì)于雷達(dá)之間的運(yùn)動(dòng)有平動(dòng)和微動(dòng)兩種情況，本文研究的彈頭目標(biāo)只考慮微動(dòng)情況下的模型，暫時(shí)忽略平動(dòng)的影響，利用EPA-BOR與電磁散射特性分析方法獲得3個(gè)目標(biāo)的回波，如圖6所示。

圖6 三個(gè)源信號(hào)Fig.6 The three source signals

隨機(jī)生成的混合矩陣為：

(14)

由此，3個(gè)源信號(hào)經(jīng)過(guò)混合矩陣得到如圖7所示的兩個(gè)混合信號(hào)。

圖7 兩個(gè)混合信號(hào)Fig.7 The two mixing signals

由于在線性瞬時(shí)混合的模型中，混合矩陣是實(shí)數(shù)矩陣，混合信號(hào)的實(shí)部包含了有效信息，因此對(duì)混合信號(hào)的實(shí)部首先進(jìn)行數(shù)據(jù)預(yù)處理再進(jìn)行密度聚類(lèi)，得到如圖8所示的聚類(lèi)結(jié)果圖。圖中一共有3個(gè)不同的簇，不在大密度區(qū)域的點(diǎn)屬于離群點(diǎn)，然后對(duì)每一個(gè)簇求其均值，得到估計(jì)的混合矩陣為：

圖8 聚類(lèi)結(jié)果Fig.8 The result of the clustering

(15)

采用文獻(xiàn)[25]中基于時(shí)頻分析的算法求得的混合矩陣為：

(16)

表2為兩種方法估計(jì)的混合矩陣的相似度和RRMSE，可見(jiàn)本文算法能夠更加準(zhǔn)確地估出混合矩陣。

表2 混合矩陣的相似度和RRMSETab.2 Similarity and RRMSE of mixing matrix

采用基于lp范數(shù)最小化的算法與采用文獻(xiàn)[26]中基于壓縮感知的算法得到的分離信號(hào)與源信號(hào)進(jìn)行對(duì)比，結(jié)果如圖9所示。

圖9 分離信號(hào)與源信號(hào)對(duì)比圖Fig.9 The contrast diagram of the separated signals and source signals

分別計(jì)算兩種算法與源信號(hào)的相似度和RRMSE，如表3所列。

表3 分離結(jié)果表Tab.3 The result table of the separation

結(jié)果表明，本文提出的算法得到的分離信號(hào)與源信號(hào)相似度較高，都達(dá)到99%以上，誤差也在一定范圍中，都在10%以下，與文獻(xiàn)[26]得到的分離結(jié)果比較，分離性能更加優(yōu)越。

4 結(jié)論

本文以空中微動(dòng)群目標(biāo)的BSS為研究目的，給出了一種基于密度聚類(lèi)和lp范數(shù)最小化相結(jié)合的高效分離方法。實(shí)驗(yàn)結(jié)果表明，本文提出的算法能夠有效應(yīng)用在彈頭目標(biāo)回波的盲分離中。利用該算法估算的混合矩陣精確度高，分離的信號(hào)與源信號(hào)相比誤差較小，在解決一些密度聚類(lèi)算法中需要事先輸入類(lèi)別個(gè)數(shù)問(wèn)題的同時(shí)，避免一些源分離最優(yōu)化算法中出現(xiàn)局部極小值的問(wèn)題，使分離性能大大提高，更加符合實(shí)際場(chǎng)景。當(dāng)群目標(biāo)之間距離較近時(shí)會(huì)存在耦合關(guān)系，因此未來(lái)的工作重點(diǎn)將解決群目標(biāo)之間的相互耦合作用，存在耦合的群目標(biāo)盲分離是一個(gè)有意義且具研究?jī)r(jià)值的問(wèn)題。本文考慮的是一種理想的場(chǎng)景，未考慮噪聲對(duì)系統(tǒng)的影響，實(shí)際中噪聲是不可忽略的，如何在有效去除噪聲的情況下達(dá)到更好的分離效果是接下來(lái)的又一重點(diǎn)工作。