恒虛警檢測信源數的方法

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

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

0 引 言

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

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

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

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

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

1 問題描述

1.1 信號模型

在譜估計[36-37]、陣列信號處理[4]、雷達成像[38]、天文學[39-40]以及其他領域[41]中,信號的觀測模型通常可以描述成如下線性模型[42]:

(1)

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

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

y=As+n

(2)

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

Ry=E[yyH]=Rs+Rn

(3)

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

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

Ry=ARsAH+σ2I

(4)

1.2 性能評價的標準

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

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

2 恒虛警的檢測法

2.1 特征值篩選

為了區分式(3)中信號與噪聲子空間中的特征值,本文首先定義一個由式(3)中相鄰特征值之差統計量構成的五維矢量如下:

(10)

由文獻[21]和文獻[30]的分析可知,理論上相鄰噪聲子空間中的特征值之差近似相等且接近于0,而信號和噪聲子空間中的特征值之差則較大。故當第i個特征值為噪聲子空間中的特征值時,序列di中的5個統計量都應該近似等于0;而當其為信號特征值時,式(10)中的5個統計量都將大于0,這樣視式(10)中不同di為不同的序列,那么就可利用K-means聚類算法將d=[d1,d2,…,dM-2]中的序列分為兩類,并分別視它們為信號和噪聲子空間。

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

步驟 1選取聚類中心。與文獻[30]的SORTE算法一樣,本文假設最小的3個特征值是屬于噪聲的。因此,將dM-2作為噪聲的聚類中心,記作μN;計算d與μN的歐氏距離,將與μN歐式距離最遠的序列作為信號的聚類中心,即

(11)

步驟 2對數據進行聚類。對每一個序列Si計算D(di,μS)和D(di,μN),如果D(di,μS)

2.2 門限計算

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

(12)

(13)

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

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

(14)

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

3 仿真實驗

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

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

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

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

表1 本文算法達到的虛警率與高斯混合模型中高斯分布個數的關系

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

表2 本文算法達到的虛警率與SNR的關系

表3 本文近似算法達到的虛警率與SNR的關系

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

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

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

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

4 結 論

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