張一迪, 王悅斌, 王培志, 楊 沁, 陸起涌, 張建秋, 李 旦

(復(fù)旦大學(xué)信息科學(xué)與工程學(xué)院, 上海 200433)

0 引 言

在生物[1]、無線通信[2]、地球物理[3]、陣列[4]、金融[5-6]等諸多信號(hào)處理任務(wù)中,預(yù)先知道待處理信號(hào)源的個(gè)數(shù)(信源數(shù))是許多信號(hào)處理算法的先驗(yàn)要求[7],如:參數(shù)化譜估計(jì)中的多信號(hào)分類(multiple signal classification,MUSIC)法[8-9]和信號(hào)參數(shù)估計(jì)的旋轉(zhuǎn)不變性技術(shù)(estimate of signal parameters by rotational invariance techniques, ESPRIT)[10-11]、盲源分離中的獨(dú)立元分析法[12]、信號(hào)降維[13]和/或特征提取中的主成分分析法[14]、參數(shù)化的時(shí)頻分析法[15]等。

近幾十年來,人們?yōu)檫@些信號(hào)處理算法發(fā)展出了一系列信源數(shù)檢測的方法,縱觀這些方法,不難發(fā)現(xiàn)它們不是基于信號(hào)觀測的特征值或特征向量,就是基于描述信號(hào)信息的信息論[16]。發(fā)展于信息論的方法主要有:赤池信息準(zhǔn)則[17-18](akaike information criterion,AIC)和最小描述長度[18-20](minimum description length,MDL)準(zhǔn)則。這兩個(gè)準(zhǔn)則都發(fā)展于信號(hào)觀測的對數(shù)似然函數(shù),以及與待估計(jì)信源數(shù)有關(guān)的懲罰項(xiàng)。分析表明:AIC不是一個(gè)一致估計(jì),即使在高信噪比(signal to noise ratio, SNR)的情況下,其結(jié)果的虛警率一直很高[21]。MDL準(zhǔn)則雖然是一個(gè)一致估計(jì),但其在低SNR或觀測數(shù)量較少的情況下,給出的高漏報(bào)率是許多應(yīng)用難以接受的[21]。此外,AIC和MDL方法對高斯白噪聲這一假設(shè),使得其在應(yīng)對非高斯噪聲時(shí),性能變得非常差[21]。針對AIC和MDL存在的這些問題,人們通過利用不同的函數(shù)來替代對數(shù)似然函數(shù)和/或改進(jìn)懲罰項(xiàng),發(fā)展出了一系列基于信息論的新方法[22-26]。盡管這些新方法在噪聲子空間中的特征值近似相等的情況下,效果較好;可是在當(dāng)觀測數(shù)據(jù)較小和/或SNR較低,以及噪聲子空間中的特征值差異較大時(shí),它們的性能都在一定程度上受到了破壞[27]?；谛盘?hào)觀測的特征向量,文獻(xiàn)[4]和文獻(xiàn)[28]報(bào)道了利用先驗(yàn)知識(shí)來檢測信源數(shù)的新方法,不過其性能一直依賴先驗(yàn)知識(shí)是否有效[27]。

最近,文獻(xiàn)[27]和文獻(xiàn)[29]報(bào)道了通過尋找信號(hào)特征值或噪聲子空間中的特征值來進(jìn)行信源數(shù)估計(jì)的方法,認(rèn)為降序排列的噪聲子空間中的特征值近似服從指數(shù)分布,如此就可利用遞歸的方法估計(jì)出一個(gè)理想噪聲子空間的特征值,當(dāng)觀測和理想噪聲子空間特征值之差且與理想噪聲子空間特征值之比大于某一閾值門限時(shí),就認(rèn)為該特征值是屬于信號(hào)的,否則就是噪聲的。該方法的主要缺點(diǎn)是:對于每一個(gè)特征值都需要采用蒙特卡羅方法計(jì)算出一個(gè)門限,文獻(xiàn)[27]報(bào)道的一個(gè)在特征值的數(shù)量為5且虛警不超過0.01時(shí)的算例表明，計(jì)算出一個(gè)門限就需要進(jìn)行16萬次蒙特卡羅仿真,且這個(gè)計(jì)算量會(huì)隨著特征值數(shù)量的增加而呈指數(shù)增加,如此大的計(jì)算量表明其應(yīng)用價(jià)值十分有限。

針對非高斯噪聲中信源數(shù)的探測問題,文獻(xiàn)[21]報(bào)道了一種對高斯白噪聲和非高斯噪聲都有效的判別函數(shù)法,它試圖利用降序排列特征值之間的斜率變化來區(qū)分信號(hào)和噪聲子空間,不過本文的研究表明:在信源數(shù)遠(yuǎn)小于觀測特征值數(shù)目時(shí),估計(jì)的虛警率居高不下?；谛盘?hào)觀測特征值的二階統(tǒng)計(jì)量(the second order statistic of eigenvalues,SORTE),文獻(xiàn)[30]給出了一種新的信源數(shù)探測法,通過計(jì)算降序排列特征值兩個(gè)相鄰特征值之差的方差來估計(jì)信源數(shù)。相較于AIC和MDL,文獻(xiàn)[30]的SORTE方法不需要估計(jì)參數(shù)來計(jì)算似然函數(shù),而與文獻(xiàn)[27]算法相比,SORTE的計(jì)算復(fù)雜度則低得多。本文將其與判別函數(shù)法進(jìn)行的比較研究發(fā)現(xiàn):雖然SORTE在信源數(shù)遠(yuǎn)小于觀測特征值數(shù)量的情況時(shí),其結(jié)果的虛警率低于判別函數(shù)法,但是當(dāng)最小的幾個(gè)特征值之差近似相等且都大于0時(shí),它會(huì)將除最小的3個(gè)特征值以外的所有特征值都判斷為信號(hào)的特征值,這意味著此時(shí)算法虛警率會(huì)高到難以滿足應(yīng)用的要求。

在目標(biāo)檢測中,許多應(yīng)用都希望檢測到假目標(biāo)的概率能保持在一個(gè)可以接受的范圍內(nèi)。比如:在雷達(dá)和聲吶的應(yīng)用中,虛警率越高意味著檢測到假目標(biāo)的數(shù)量可能越多,而恒虛警檢測則是指其中假目標(biāo)所占的比例是固定的[31]?？墒?文獻(xiàn)中報(bào)道的信源數(shù)檢測方法都不是恒虛警的,因此存在產(chǎn)生過高虛警率的可能性,這也就意味著其在滿足雷達(dá)和/或聲納目標(biāo)檢測的應(yīng)用要求時(shí)存在問題。針對存在的這一問題,本文提出了一種恒虛警檢測信源數(shù)的方法,該方法首先將一個(gè)M×M階觀測協(xié)方差矩陣相鄰特征值之差的統(tǒng)計(jì)量,定義成一個(gè)五維矢量序列,這樣這個(gè)序列就有M-2個(gè),利用K均值(K-means)聚類算法[32-35]將這M-2個(gè)五維矢量序列分成兩類,并分別視其為信號(hào)和噪聲兩個(gè)子空間。當(dāng)將噪聲子空間中的特征值序列描述成一個(gè)統(tǒng)計(jì)分布,并通過期望最大(expectation maximization,EM)算法估計(jì)出這個(gè)統(tǒng)計(jì)分布時(shí),奈曼-皮爾遜(Neyman-Pearson,NP)假設(shè)檢驗(yàn)就可利用這個(gè)分布來對信源數(shù)進(jìn)行恒虛警檢測?？紤]到EM算法估計(jì)統(tǒng)計(jì)分布和蒙特卡羅算法計(jì)算門限時(shí)計(jì)算復(fù)雜度過高,本文也給出了一個(gè)不管噪聲子空間中的特征值服從何分布,且可大大降低計(jì)算復(fù)雜度的近似NP假設(shè)檢驗(yàn)方法,這樣就進(jìn)一步簡化了本文方法的計(jì)算并增強(qiáng)了它的實(shí)用性。數(shù)值仿真的結(jié)果在驗(yàn)證提出方法的有效性的同時(shí),也優(yōu)于其他方法。

1 問題描述

1.1 信號(hào)模型

在譜估計(jì)[36-37]、陣列信號(hào)處理[4]、雷達(dá)成像[38]、天文學(xué)[39-40]以及其他領(lǐng)域[41]中,信號(hào)的觀測模型通?？梢悦枋龀扇缦戮€性模型[42]:

(1)

式中：y=[y1,y2,…yM]表示觀測信號(hào)；ωk和sk為第k個(gè)信號(hào)分量的未知參數(shù)；a(·)是一個(gè)具有未知參數(shù)ωk的已知函數(shù)；K是未知的信源數(shù)；n為噪聲。

重寫式(1)為如下矩陣形式:

y=As+n

(2)

式中：A=[a(ω1),a(ω2),…,a(ωk)]；s=[s1(t),s2(t),…sk(t)]T。此時(shí),式(2)的協(xié)方差矩陣可以表示為[27]

Ry=E[yyH]=Rs+Rn

(3)

式中：E[·]表示期望運(yùn)算。在無噪聲的理想情況下,有Ry=Rs,也就是說信源數(shù)是觀測信號(hào)y協(xié)方差矩陣Ry的秩,即K=rank(Ry)。而在應(yīng)用中,由于噪聲的存在,矩陣Ry一般是滿秩的,因此信源數(shù)K≤rank(Ry),這意味著一般不能通過觀測信號(hào)協(xié)方差矩陣的秩來確定信源數(shù)[27]。

在理想白噪聲情況下,式(3)可以改寫為[18,27]

Ry=ARsAH+σ2I

(4)

1.2 性能評價(jià)的標(biāo)準(zhǔn)

本文將用虛警率、漏報(bào)率、檢測能力以及錯(cuò)誤率來評價(jià)算法的性能。用K表示實(shí)際信源數(shù),K*表示估計(jì)信源數(shù),M表示觀測數(shù),這樣虛警率、漏報(bào)率、檢測能力以及錯(cuò)誤率的定義以及計(jì)算方法就可以分別給出如下[35]。

檢測能力(PD)則表示信號(hào)特征值判別為信號(hào)的概率,計(jì)算公式為PD=1-PM。

2 恒虛警的檢測法

2.1 特征值篩選

為了區(qū)分式(3)中信號(hào)與噪聲子空間中的特征值,本文首先定義一個(gè)由式(3)中相鄰特征值之差統(tǒng)計(jì)量構(gòu)成的五維矢量如下:

(10)

由文獻(xiàn)[21]和文獻(xiàn)[30]的分析可知,理論上相鄰噪聲子空間中的特征值之差近似相等且接近于0,而信號(hào)和噪聲子空間中的特征值之差則較大。故當(dāng)?shù)趇個(gè)特征值為噪聲子空間中的特征值時(shí),序列di中的5個(gè)統(tǒng)計(jì)量都應(yīng)該近似等于0;而當(dāng)其為信號(hào)特征值時(shí),式(10)中的5個(gè)統(tǒng)計(jì)量都將大于0,這樣視式(10)中不同di為不同的序列,那么就可利用K-means聚類算法將d=[d1,d2,…,dM-2]中的序列分為兩類,并分別視它們?yōu)樾盘?hào)和噪聲子空間。

令CN和CS分別表示噪聲和信號(hào)子空間,μN(yùn)和μs分別表示其所對應(yīng)的聚類中心,則由如下K-means算法就可分別得到式(3)的噪聲和信號(hào)子空間:

步驟 1選取聚類中心。與文獻(xiàn)[30]的SORTE算法一樣,本文假設(shè)最小的3個(gè)特征值是屬于噪聲的。因此,將dM-2作為噪聲的聚類中心,記作μN(yùn);計(jì)算d與μN(yùn)的歐氏距離,將與μN(yùn)歐式距離最遠(yuǎn)的序列作為信號(hào)的聚類中心,即

(11)

步驟 2對數(shù)據(jù)進(jìn)行聚類。對每一個(gè)序列Si計(jì)算D(di,μS)和D(di,μN(yùn)),如果D(di,μS)

2.2 門限計(jì)算

為了在給定虛警率的條件下檢測信源數(shù),就需要根據(jù)NP準(zhǔn)則來設(shè)計(jì)判決門限δn,這樣特征值中滿足不等式λi≥δn的數(shù)量,就可視為給定虛警率的信源數(shù)[35]。通過NP準(zhǔn)則來計(jì)算判別的門限,則需要知道噪聲子空間中的特征值所服從分布的概率密度函數(shù),而概率密度函數(shù)的高斯混合模型則可描述任一概率分布函數(shù)[43],因此假定式(3)噪聲子空間中的特征值的概率密度函數(shù)可由如下高斯混合模型描述:

(12)

(13)

考慮到使用EM算法估計(jì)高斯混合模型的參數(shù)以及用蒙特卡羅方法計(jì)算門限時(shí)計(jì)算復(fù)雜度較高,因此,假設(shè)噪聲子空間中的特征值是概率分布,都可由高斯分布近似。此時(shí),近似的NP假設(shè)檢驗(yàn)的判別門限為[35]

δn=μN(yùn)+Q-1(1-α)σN

(14)

式中:μN(yùn)和σN分別表示噪聲子空間中特征值的均值和方差,Q-1表示正態(tài)分布累積分布函數(shù)的反函數(shù)。這樣,就只需要計(jì)算噪聲子空間中特征值的均值和方差就可以得到判別的門限。雖然這樣處理大大減小了本文算法的計(jì)算復(fù)雜度,但是當(dāng)噪聲子空間中的特征值的概率分布偏離假設(shè)的高斯分布較大時(shí),虛警率的誤差會(huì)變大。即使如此,下面的仿真研究也表明:其性能依然優(yōu)于文獻(xiàn)最新報(bào)道的方法。

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

為了驗(yàn)證本文算法的有效性,本節(jié)進(jìn)行了7組仿真實(shí)驗(yàn),分別研究了提出算法控制虛警率的能力,檢測能力、錯(cuò)誤率和虛警率,及其與特征值的長度、SNR以及目標(biāo)數(shù)量之間的關(guān)系,并將提出方法與文獻(xiàn)中最新報(bào)道的算法進(jìn)行了對比。為了后面描述方便,在不混淆的情況下,將用高斯混合模型描述噪聲子空間中特征值的統(tǒng)計(jì)分布得到的算法作為本文算法;而本文近似算法則代表假設(shè)噪聲子空間中特征值的統(tǒng)計(jì)分布是高斯分布,并用式(13)計(jì)算門限的算法。

第1組仿真實(shí)驗(yàn),將研究本文算法控制虛警率的誤差與信源數(shù),以及噪聲子空間中的特征值數(shù)N的關(guān)系,這里噪聲子空間中的特征值數(shù)目與信號(hào)特征值數(shù)目之和簡稱為觀測數(shù)。在仿真中,對高斯混合分布的混合個(gè)數(shù),及其均值和方差的隨機(jī)變化進(jìn)行了1萬次蒙特卡羅仿真,對其結(jié)果求均值得到的虛警率最大誤差如圖1所示。從圖1中可以看出在噪聲子空間中的特征值數(shù)量N≥30時(shí),虛警率最大誤差幾乎不隨信源數(shù)而改變。

第2組仿真實(shí)驗(yàn)將研究本文算法控制虛警率的誤差與信源數(shù)以及SNR的關(guān)系。仿真中,在噪聲子空間中的特征值數(shù)量為30的情況下對高斯混合分布的混合個(gè)數(shù),以及它們均值和方差的隨機(jī)變化的噪聲進(jìn)行了1萬次蒙特卡羅仿真,對其結(jié)果求均值后得到的虛警率最大誤差如圖2所示。圖2(a)的仿真實(shí)驗(yàn)結(jié)果表明,本文算法在SNR大于0 dB時(shí),虛警率最大誤差幾乎不隨信源數(shù)而改變。而圖2(b)的仿真實(shí)驗(yàn)結(jié)果則表明,本章近似算法在SNR大于0 dB且不超過5 dB時(shí),虛警率最大誤差雖然隨信源數(shù)的變化而有所改變,但其改變的數(shù)值較小。當(dāng)SNR大于5 dB時(shí),虛警率誤差隨SNR的增大而增加。這是因?yàn)殡S著SNR增大,噪聲子空間中的特征值近似相等,此時(shí)其實(shí)際的統(tǒng)計(jì)分布與高斯分布相差較大。

第3組仿真實(shí)驗(yàn)研究了本文算法給定的虛警率與實(shí)際可達(dá)到的虛警率之間的誤差,同時(shí)研究了高斯混合模型中高斯分布的個(gè)數(shù)c對其誤差的影響。在仿真中,噪聲子空間中的特征值數(shù)量設(shè)為30,在信源數(shù)隨機(jī)變化的情況下,對高斯混合分布的混合個(gè)數(shù),及其均值和方差的隨機(jī)變化的噪聲進(jìn)行了1萬次蒙特卡羅仿真,對其結(jié)果求均值后得到了表1的數(shù)據(jù)。從表1的數(shù)據(jù)中可以看出,本文算法控制虛警率的誤差在±5%,且隨著高斯混合模型中高斯分布數(shù)量的增加,虛警率誤差呈減小趨勢,當(dāng)高斯分布個(gè)數(shù)大于等于5時(shí),高斯分布個(gè)數(shù)對虛警率誤差的影響幾乎為0。

表1 本文算法達(dá)到的虛警率與高斯混合模型中高斯分布個(gè)數(shù)的關(guān)系

第4組仿真實(shí)驗(yàn)研究了本文算法給定的虛警率與實(shí)際可達(dá)到的虛警率之間的誤差,同時(shí)研究了SNR對其誤差的影響。在仿真中,噪聲子空間中的特征值數(shù)量設(shè)為30,在信源數(shù)隨機(jī)變化的情況下,對高斯混合分布的混合個(gè)數(shù),及其均值和方差的隨機(jī)變化的噪聲進(jìn)行了1萬次蒙特卡羅仿真,對其結(jié)果求均值后得到了表2和表3的數(shù)據(jù)。從表2的數(shù)據(jù)中可以看出,當(dāng)本文算法用高斯混合模型描述噪聲子空間中的特征值的統(tǒng)計(jì)分布時(shí),算法控制虛警率的誤差在±5%,且SNR對虛警率的誤差的影響可以忽略。而由表3的數(shù)據(jù)可知,當(dāng)本文算法假設(shè)噪聲子空間中的特征值服從高斯分布時(shí),算法控制虛警率的誤差會(huì)隨著SNR的增加而變大。這是因?yàn)殡S著SNR的增加,噪聲子空間中的特征值近似相等,此時(shí)其實(shí)際的統(tǒng)計(jì)分布與高斯分布相差較大。

表2 本文算法達(dá)到的虛警率與SNR的關(guān)系

表3 本文近似算法達(dá)到的虛警率與SNR的關(guān)系

第5組仿真實(shí)驗(yàn)將對判別函數(shù)法、SORTE算法、本文算法和本文近似算法,在不同信源數(shù)的情況下進(jìn)行對比研究。在仿真中,觀測數(shù)M=50,SNR為0 dB,虛警率給定為α=0.05。圖3給出了這3種算法的虛警率、錯(cuò)誤率和檢測能力分別在10 000次蒙特卡羅仿真中的平均結(jié)果。

圖3(a)仿真結(jié)果表明:SORTE算法的虛警率隨著信源數(shù)的增加而有增大的趨勢;判別函數(shù)法的虛警率隨著信源數(shù)的增加而減小;本文算法與本文近似算法的虛警率幾乎保持恒定,與給定的虛警率的誤差不超過±5%。圖3(b)表明:當(dāng)信源數(shù)小于8時(shí),本文算法及本文近似算法的錯(cuò)誤率與SORTE算法相差不大,小于判別函數(shù)法;而當(dāng)信源數(shù)大于8時(shí),本文算法及本文近似算法的錯(cuò)誤率小于SORTE算法,略高于判別函數(shù)法。但是在只有一個(gè)目標(biāo)時(shí),本文算法及本文近似算法的錯(cuò)誤率遠(yuǎn)小于判別函數(shù)法。圖3(c)當(dāng)信源數(shù)小于15時(shí),本文算法及本文近似算法的檢測能力略差于判別函數(shù)法,好于SORTE算法;而當(dāng)信源數(shù)大于15時(shí),本文算法與本文近似算法的檢測能力略好于判別函數(shù)法和SORTE算法。

第6組仿真實(shí)驗(yàn)將研究判別函數(shù)法、SORTE算法和本文算法與SNR的關(guān)系。在仿真中,觀測數(shù)為35,信源數(shù)為3,其中包含2個(gè)大目標(biāo)1個(gè)小目標(biāo),且小目標(biāo)的幅度為大目標(biāo)的-5 dB,虛警率給定為α=0.05。圖4給出了這3種算法的虛警率、錯(cuò)誤率和檢測能力分別在10 000次蒙特卡羅仿真中的平均結(jié)果。圖4(a)仿真結(jié)果表明:SNR對SORTE算法的虛警率影響幾乎可以忽略;判別函數(shù)法的虛警率則隨SNR的增加而減小;而本文算法與其近似算法的虛警率幾乎保持恒定,與給定的虛警率的誤差不超過±5%。圖4(b)則表明:本文算法及本文近似算法的錯(cuò)誤率小于判別函數(shù)法和SORTE算法。圖4(c)則表明:本文算法及本文近似算法的檢測能力好于SORTE算法,與判別函數(shù)法幾乎一樣。但由圖4(b)可知,判別函數(shù)法的錯(cuò)誤率是最高的;這是由于該算法在目標(biāo)數(shù)量較少時(shí)虛警率過高造成的。

第7組仿真實(shí)驗(yàn)將研究判別函數(shù)法、SORTE算法和本文算法與觀測數(shù)的關(guān)系。在仿真中,SNR為0 dB,信源數(shù)為5,其中包含3個(gè)大目標(biāo)和2個(gè)小目標(biāo),且小目標(biāo)的幅度分別為大目標(biāo)的-5 dB和-10 dB,虛警率給定為α=0.05。圖5給出了這3種算法的虛警率、錯(cuò)誤率和檢測能力分別在1萬次蒙特卡羅仿真中的平均結(jié)果。圖5(a)仿真結(jié)果表明:SORTE算法的虛警率隨觀測數(shù)增加而減小;判別函數(shù)法的虛警率隨著觀測數(shù)的增加而增加;本文算法與其近似算法的虛警率幾乎保持恒定,與給定的虛警率的誤差不超過±5%。同時(shí),可以看出隨著觀測數(shù)增加,本文算法及本文近似算法控制虛警率的誤差變小。圖5(b)表明:本文算法及本文近似算法的錯(cuò)誤率小于判別函數(shù)法和SORTE算法。圖5(c)則表明:當(dāng)目標(biāo)數(shù)量較少且含有小目標(biāo)時(shí),判別函數(shù)法的檢測能力最好,但由圖5(b)可知,此時(shí)該種算法的錯(cuò)誤率是最高的,這是由于其虛警率較高而導(dǎo)致的;SORTE算法則只能檢測出3個(gè)大目標(biāo)。而本文算法及本文近似算法的檢測能力則好于SORTE算法,略低于判別函數(shù)法。需要強(qiáng)調(diào)的是,本文算法及本文近似算法的虛警率是最低的,且控制虛警率誤差的能力最好,這就意味著其可以通過提高虛警率來提高本文算法的檢測能力。

4 結(jié) 論

針對目前已有的信源數(shù)檢測算法不能控制虛警率的問題,本文提出了一種恒虛警檢測信源數(shù)的方法。該方法利用本文定義的觀測信號(hào)協(xié)方差矩陣相鄰特征值之差統(tǒng)計(jì)量所構(gòu)成的五維矢量序列,可利用K-means聚類算法將其分成兩類。當(dāng)分別視其為信號(hào)和噪聲兩個(gè)子空間,且將噪聲子空間對應(yīng)的特征值序列描述成一個(gè)統(tǒng)計(jì)分布,并通過EM算法估計(jì)出這個(gè)統(tǒng)計(jì)分布時(shí),NP假設(shè)檢驗(yàn)就可利用這個(gè)分布來對信源數(shù)進(jìn)行恒虛警檢測。為了提高算法效率和實(shí)用性,本文也給出了一種可大大降低計(jì)算復(fù)雜度的近似奈曼-皮爾遜假設(shè)檢驗(yàn)方法。數(shù)值仿真的結(jié)果在驗(yàn)證提出方法的有效性的同時(shí),也表明其性能優(yōu)于判別函數(shù)法和SORTE算法。