非均勻噪聲下基于BOOTSTRAP和特征空間投影的信源數(shù)估計(jì)

郭立民，馮 凱

（哈爾濱工程大學(xué) 信息與通信工程學(xué)院，哈爾濱150001）

0 引 言

陣列信號(hào)處理在許多應(yīng)用領(lǐng)域有著重要作用。信源數(shù)檢測(cè)作為信號(hào)參數(shù)估計(jì)的先決條件，是陣列信號(hào)處理中至關(guān)重要的問(wèn)題［1－2］。最早的信源數(shù)檢測(cè)問(wèn)題是采用假設(shè)檢驗(yàn)，需要人為地設(shè)置門限［3］，因而人們提出了預(yù)先設(shè)定包含門限置信區(qū)間的檢測(cè)方法［4］。經(jīng)典的信源數(shù)檢測(cè)方法以AIC和MDL 信息論準(zhǔn)則為基礎(chǔ)［5］，AIC 準(zhǔn)則缺乏一致性，MDL 準(zhǔn)則是一致估計(jì)，但在低信噪比和小快拍數(shù)的情況下不能保證檢測(cè)性能。在文獻(xiàn)［6，7］中對(duì)MDL 準(zhǔn)則做出一些修改來(lái)提高其檢測(cè)性能。近年來(lái)，為了適應(yīng)更加真實(shí)的信號(hào)環(huán)境，提出了許多新方法：文獻(xiàn)［8］使用了觀測(cè)值的概率密度函數(shù)，在信息論準(zhǔn)則的框架下從樣本協(xié)方差矩陣獲得樣本特征值，其缺點(diǎn)是在實(shí)際陣列系統(tǒng)中數(shù)據(jù)分布往往是未知的。一般的信源數(shù)檢測(cè)方法都假設(shè)噪聲為空間均勻的，噪聲的均勻性能夠完全利用數(shù)據(jù)協(xié)方差矩陣的特征值所包含的信息，然后實(shí)際的陣列噪聲往往是空間非均勻的，即各傳感器信道噪聲不一致時(shí)，使得大部分的算法失效。在文獻(xiàn)［9］中提到了空間非均勻噪聲下信源數(shù)檢測(cè)問(wèn)題的模型。蓋氏圓方法［10－12］能夠適用于空間非均勻噪聲的環(huán)境，但其在信噪比較低時(shí)常常會(huì)產(chǎn)生欠估計(jì)或過(guò)估計(jì)，且需要對(duì)（M－1）×（M－1）階的矩陣特征值分解，而在波達(dá)方向估計(jì)中需要M×M 階特征值分解，因此在獲得波達(dá)方向是需要做兩次矩陣特征值分解，運(yùn)算量過(guò)大。

在許多實(shí)際的系統(tǒng)中，高斯數(shù)據(jù)和大樣本數(shù)量的假設(shè)是無(wú)法滿足的［13］，所以本文提出了一種基于假設(shè)檢驗(yàn)的替代方法。提出的假設(shè)檢驗(yàn)以特征空間投影為基礎(chǔ)。零假設(shè)下的檢驗(yàn)統(tǒng)計(jì)量的分布是使 用Bootstrap 方 法 來(lái) 估 計(jì) 的［14－15］，這 種 方法對(duì)高斯分布和大部分非高斯分布都是有效的，且其優(yōu)點(diǎn)是估計(jì)在確定樣本數(shù)量下的經(jīng)驗(yàn)分布，而不是漸進(jìn)分布。因此，提出新方法的優(yōu)點(diǎn)是可以解決不滿足高斯分布和漸進(jìn)分布條件下的問(wèn)題。仿真結(jié)果表明，所提出的方法在少的樣本數(shù)量和低信噪比下，有較好的估計(jì)性能。

1 數(shù)據(jù)模型

設(shè)N 個(gè)遠(yuǎn)場(chǎng)窄帶信號(hào)入射到由M 個(gè)陣元組成的均勻線陣，入射角為θ1，…，θN，陣元間距為d＝λ／2。信源數(shù)N 為未知估計(jì)量，則線陣接收的陣列信號(hào)為［9］：

式中：X（）t，S（）t 和N（）t 為均值為零的隨機(jī)過(guò)程，為陣列導(dǎo)向矩陣，S（）t 是信號(hào)矢量，N（）t 是噪聲矢量。

加性噪聲假設(shè)為零均值的空間非均勻白噪聲，其協(xié)方差矩陣滿足下式：

這里（）·H為矩陣的共軛轉(zhuǎn)置算子，空間非均勻噪聲各傳感器能量各不相同信號(hào)和噪聲假設(shè)是不相關(guān)的，且各自分布未知。陣列協(xié)方差矩陣為：

2 信源數(shù)檢驗(yàn)

2.1 非均勻噪聲下的MDL準(zhǔn)則

由于傳感器噪聲能量各不相同，令陣列的任意元素失效來(lái)減輕各傳感器噪聲的作用。簡(jiǎn)單而不失一般性，第M 個(gè)傳感器失效，協(xié)方差矩陣進(jìn)行分塊：

取M－1維方差的特征空間（即特征矩陣E，滿足EEΗ＝I 且RM＝EΣ′EΗ）構(gòu)成酉變化矩陣［11］：

這樣酉變化之后的矩陣為：

矢量c的元素cm幅度值是蓋氏圓半徑［11］，其指出了特征值多樣性和特征矢量對(duì)應(yīng)空間［16］：

從上式可以清楚地看出，兩個(gè)不同的圓盤，前N個(gè)半徑ρ1…ρN代表了信號(hào)子空間，最小的M－1－N 半徑代表了噪聲子空間。

當(dāng)噪聲為空間非均勻的情況下，不能通過(guò)協(xié)方差矩陣的特征值排序來(lái)進(jìn)行信號(hào)源檢測(cè)，我們用所包含的信息來(lái)代替。容易得到觀測(cè)數(shù)據(jù)的負(fù)對(duì)數(shù)似然方程為［17］：

式（5）定義的酉矩陣TTΗ＝I，我們可以用協(xié)方差變換矩陣RT代替式（8）中的R，省略獨(dú)立于N 的部分，最終我們推導(dǎo)得出的負(fù)對(duì)數(shù)似然函數(shù)變?yōu)椋?/p>

根據(jù)文獻(xiàn)［5］，罰函數(shù)可以用來(lái)進(jìn)行信息論準(zhǔn)則的信源數(shù)檢測(cè)。一般，MDL 準(zhǔn)則可以表達(dá)如下：

MDL準(zhǔn)則的罰函數(shù)被定義為：

2.2 基于蓋世圓半徑的假設(shè)檢驗(yàn)

在實(shí)際中，如果所有的蓋氏圓半徑都對(duì)應(yīng)于噪聲空間，統(tǒng)計(jì)量T1i（半徑平方累計(jì)和）和T2i（算術(shù)均值和幾何均值的差異）將非常?。ń咏?），否則會(huì)比零值大很多。根據(jù)這種變化，進(jìn)行以下一系列獨(dú)立的假設(shè)檢驗(yàn)［19］：

檢測(cè)從H0開(kāi)始，當(dāng)接收H0時(shí)意味著沒(méi)有信號(hào)源存在，信號(hào)源數(shù)目為N ＝0。如果H0被拒絕，那么相對(duì)的）被接受，意味著至少存在一個(gè)信號(hào)，然而沒(méi)有信息能夠說(shuō)明具體存在幾個(gè)信號(hào)。因此，通過(guò)逐步進(jìn)行假設(shè)Ηi，最大的蓋氏圓半徑的作用被減小。當(dāng)某一假設(shè)被接受或者到達(dá)HM－2時(shí)，檢驗(yàn)停止。

對(duì)于給定的顯著水平α，其含義為當(dāng)前沒(méi)有信號(hào)存在時(shí)，正確檢測(cè)數(shù)信號(hào)源數(shù)目N ＝0的概率被維持在1－α。因此，每個(gè)假設(shè)都在顯著水平α下被檢驗(yàn)。

2.3 提出的新方法

我們知道R 的M－N 個(gè)最小的特征值是噪聲方差的估計(jì)值，為了減小噪聲的貢獻(xiàn)，執(zhí)行降噪運(yùn)算，取RX＝R－γ2I，其中γ2為陣列信號(hào)樣本協(xié)方差矩陣最小的兩個(gè)特征值的平均值，即γ2＝

由式（3）得出：

式中：bi是A 的行矢量，陣列協(xié)方差矩陣可表示為：

式中：Σ和UC分別為陣列信號(hào)協(xié)方差矩陣特征值組成的對(duì)角矩陣和特征矢量矩陣。可以證明，經(jīng)過(guò)降噪處理的陣列信號(hào)協(xié)方差矩陣為：

式中：Σ′為陣列信號(hào)的特征值減去噪聲功率估計(jì)γ2后組成的對(duì)角陣；UC與降噪前相同，UC可以表示為，其中ui是RX的特征矢量，RX的特征值從大到小排列，ui按與其對(duì)應(yīng)的特征值的大小排列，求RX在UC上的投影UΗCRX，其i列為：

由于信號(hào)特征子空間與導(dǎo)向矢量在同一空間內(nèi)，

各信號(hào)的功率常常是不等的，為了提高估計(jì)性能，對(duì)信號(hào)功率進(jìn)行歸一化以消除各個(gè)信號(hào)功率不同對(duì)信源數(shù)估計(jì)產(chǎn)生的影響。投影還可以表示為：

由上式，功率歸一化后，RX在UC上的投影就是，因此我們可以用UC估計(jì)目標(biāo)數(shù)。ui表示為：第m 行元素模的均值為：

如果特征矢量對(duì)應(yīng)于噪聲空間，則檢測(cè)統(tǒng)計(jì)量T1i和T2i接近于零。進(jìn)行下面的一系列獨(dú)立的假設(shè)檢驗(yàn)。實(shí)際估計(jì)信源數(shù)從Η0開(kāi)始，當(dāng)前假設(shè)被拒絕時(shí)進(jìn)行下一個(gè)假設(shè)，直到接受為止。

這里顯著水平與式（14）的假設(shè)檢驗(yàn)的顯著水平α相同，表示當(dāng)前沒(méi)有信號(hào)存在時(shí)，正確檢驗(yàn)出N＝0的概率被維持在1－α。

2.4 Bootstrap對(duì)于分布的估計(jì)

進(jìn)行假設(shè)檢驗(yàn)的顯著度估計(jì)需要檢驗(yàn)統(tǒng)計(jì)量的空域分布知識(shí)，Bootstrap方法用來(lái)進(jìn)行零假設(shè)下的非參數(shù)估計(jì)［18］。使用Bootstrap方法可以避免對(duì)數(shù)據(jù)分布做出假設(shè)。如前文提到的，這個(gè)優(yōu)點(diǎn)是非常重要的。Bootstrap方法的原理是數(shù)據(jù)樣本代表了真實(shí)分布的經(jīng)驗(yàn)估計(jì)。因此用Bootstrap 方法重新采樣的數(shù)據(jù)集進(jìn)行推理。Bootstrap方法避免了小樣本情況下的漸進(jìn)估計(jì)的錯(cuò)誤。類似地，文獻(xiàn)［20］強(qiáng)調(diào)了改善復(fù)雜的檢驗(yàn)統(tǒng)計(jì)量的經(jīng)驗(yàn)分布的估計(jì)問(wèn)題，例如對(duì)更少的樣本進(jìn)行重采樣。

這里Ti為從數(shù)據(jù)X（）t 得到的檢驗(yàn)統(tǒng)計(jì)量，為從重采樣樣本得到的檢驗(yàn)統(tǒng)計(jì)量。

文獻(xiàn)［20，21］給出假設(shè)檢驗(yàn)的顯著度定義如下：

這里I（）· 為指標(biāo)函數(shù)。因此，檢驗(yàn)量為α，從i＝0開(kāi)始，如果Pi≥α，那么Ηi被接受，否則，設(shè)i＝i＋1并且重復(fù)檢驗(yàn)。

3 仿真結(jié)果

線性均勻陣列陣元數(shù)目為M ＝5，真實(shí)信源數(shù)為N ＝2，快拍數(shù)為L(zhǎng)＝100，重采樣次數(shù)B＝100，檢驗(yàn)量α＝2%信號(hào)波達(dá)方向?yàn)棣?＝，信噪比SNR＝10dB，且最差噪聲比WNPR＝10。

這里SNR 和WNPR 定義如下：

蒙特卡洛實(shí)驗(yàn)次數(shù)為1000，除了特殊說(shuō)明，假設(shè)信號(hào)和噪聲為高斯分布。

圖1為檢測(cè)器性能與信噪比的關(guān)系對(duì)比圖。信噪比從－15dB變化到20dB，快拍數(shù)為100，可以看出，檢測(cè)器性能隨著信噪比的增大而提高。在非均勻噪聲的情況下，AIC 和MDL 信息論準(zhǔn)則失效。在低信噪比下，所提方法的性能優(yōu)于GBootstrap算法和NU－MDL 算法，同時(shí)在高信噪比下，新方法的檢測(cè)器性能略好?；贐ootstrap方法的檢測(cè)器因不需要任何數(shù)據(jù)分布的先驗(yàn)知識(shí)而表現(xiàn)優(yōu)異。

圖2為拉普拉斯信號(hào)和高斯噪聲，圖3為高斯信號(hào)和拉普拉斯噪聲，可以看出，當(dāng)信號(hào)或者噪聲不滿足高斯分布時(shí)，NU－MDL的性能并沒(méi)有在高斯分布的情況下表現(xiàn)出來(lái)?；贐ootstrap方法的檢測(cè)器不需要對(duì)數(shù)據(jù)分布做出假設(shè)，當(dāng)信噪比增加時(shí)，性能優(yōu)于NU－MDL。新算法因?qū)﹃嚵袇f(xié)方差矩陣降噪處理，應(yīng)用信號(hào)子空間和噪聲子空間的正交性，易于表征投影的大小區(qū)分信號(hào)和噪聲的貢獻(xiàn)，所以性能優(yōu)于GBootstrap算法。新算法的檢測(cè)性能同高斯分布的情況下基本相同，說(shuō)明了檢測(cè)器的魯棒性較好。

圖1 在非均勻噪聲下改變信噪比的檢測(cè)器性能Fig.1 Probability of detection against SNR

圖2 在拉普拉斯信號(hào)下改變信噪比的檢測(cè)器性能Fig.2 Probability of detection with Laplacian signals

圖3 在拉普拉斯噪聲下改變信噪比的檢測(cè)器性能Fig.3 Probability of detection with Laplacian noise

圖4 為信源功率不相等的情況下，檢測(cè)概率隨信噪比的變化。兩個(gè)信號(hào)的功率比為1∶2，信噪比以最小信號(hào)功率為基準(zhǔn)，由圖可見(jiàn)，在信源功率不相等的情況下，新方法依然好于其他兩種方法。

圖5說(shuō)明了快拍數(shù)對(duì)檢測(cè)器性能的作用。與預(yù)期的一樣，在快拍數(shù)變化時(shí)，基于Bootstrap方法的檢測(cè)器性能隨著快拍數(shù)增加而提高，當(dāng)快拍數(shù)達(dá)到60時(shí)，正確檢測(cè)概率不再變化，趨于穩(wěn)定。

圖4 信源功率不等情況下改變信噪比的檢測(cè)器性能Fig.4 Probability of detection with signal’s power different

圖5 改變快拍數(shù)的檢測(cè)器性能Fig.5 Probability of detection against sample size

圖6 改變WNPR 時(shí)的檢測(cè)器性能Fig.6 Probability of detection against WNPR

圖6 說(shuō)明當(dāng)WNPR 變化時(shí)的檢測(cè)性能。這里假設(shè)最小噪聲功率為0.1，最大的噪聲功率為0.1×WNPR，其他噪聲功率在最小噪聲功率和最大噪聲功率之間隨機(jī)取值。令WNPR 的值從1變化到45。注意到WNPR 接近1，此時(shí)最大噪聲功率和最小噪聲功率相等，相當(dāng)于噪聲為空間均勻的，所以檢測(cè)概率最高。從圖中可以看出NU－MDL和GBootstrap方法檢測(cè)器T1對(duì)非均勻噪聲的敏感性。這里為高斯數(shù)據(jù)，所提出的新檢驗(yàn)方法性能優(yōu)于NU－MDL 算法和GBootstrap算法，當(dāng)WNPR 小于15時(shí)，檢測(cè)器T1的檢測(cè)概率接近于100%，當(dāng)WNPR 大于20時(shí)，新方法的正確檢測(cè)概率不會(huì)隨著WNPR 的變化而變化。

圖7說(shuō)明了信號(hào)角度分辨率對(duì)檢測(cè)性能的影響。這里令θ1＝10°，并且θ2設(shè)置為從10°逐漸變化到26°。檢測(cè)性能隨著角度分辨率的增加而提高，注意到，角度分辨率大于6°時(shí)，用所提出的方法的檢測(cè)器T1是一致估計(jì)。

圖7 改變角度分辨率的檢測(cè)器性能Fig.7 Probability of detection against the angular resolution

圖8 說(shuō)明陣元數(shù)目M 和快拍數(shù)L 對(duì)檢測(cè)器T1性能的聯(lián)合作用?？炫臄?shù)從20變化到120，通過(guò)不同的陣元數(shù)目來(lái)衡量其性能。一般來(lái)說(shuō)，檢測(cè)器的性能隨著L 增加而得到改善。另一個(gè)重要的觀測(cè)結(jié)果是M 的作用，實(shí)際上，隨著X（t）維度的增加，Bootstrap需要更大的“最小”快拍數(shù)來(lái)估計(jì)不同變量的經(jīng)驗(yàn)分布。因此，在L ≤70，M＝8的檢測(cè)概率要小于M ＝6。另一方面，在某一個(gè)數(shù)值（比如L ＝80），Bootstrap可給出更好的結(jié)果，當(dāng)增加M 和L，其聯(lián)合作用明顯地改善了檢測(cè)器的檢測(cè)概率。

圖8 陣元和快拍數(shù)的關(guān)系以及檢測(cè)器T1 的性能Fig.8 Relation between M and Land effect on the T1detector’s performance.

4 結(jié)束語(yǔ)

文章提出了一種在空間非均勻噪聲下的信號(hào)檢測(cè)模型。這種模型基于信號(hào)空間和噪聲空間的正交性，提出了一系列假設(shè)檢驗(yàn)，并且數(shù)據(jù)分布的先驗(yàn)知識(shí)未知。所提出的檢驗(yàn)統(tǒng)計(jì)量是基于特征空間的投影，其在空域下的分布用Bootstrap 方法進(jìn)行估計(jì)。仿真結(jié)果表明，所提出的方法在少量樣本的情況下，性能優(yōu)于其他算法，尤其當(dāng)數(shù)據(jù)為非高斯數(shù)據(jù)，或者在低信噪比等不理想的條件下。當(dāng)重采樣次數(shù)B 增加時(shí)，計(jì)算量增大，然而仍在計(jì)算機(jī)可處理的能力范圍內(nèi)。因此，提出的檢測(cè)模型在信源數(shù)估計(jì)問(wèn)題中，性能優(yōu)于其他算法。

［1］Zoubir A M.Bootstrap methods for model selection［J］.International Journal of Electronics and Communications，1999，53（6）：386－392.

［2］Zoubir A M，Iskander D R.Bootstrap Techniques for Signal Processing［M］.Cambridge：Cambridge University Press，2004.

［3］Lawley D.Tests of significance for the latent roots of covariance and correlation matrices［J］.Biometrika，1956，43：128－136.

［4］Chen W，Wong K M，Reilly J P.Detection of the number of signals：apredicted eigen－threshold approach［J］.IEEE Trans on Signal Processing，1991，39（5）：1088－1098.

［5］Wax M，Kailath T.Detection of signals by information theoretic criteria［J］.IEEE Transactions Acoustics Speech Signal Processing，1985，33（2）：387－392.

［6］Fishler E，Messer H.Order statistics approach for determining the number of sources using an array of sensors［J］.IEEE Signal Processing Letter，1999，6：179－182.

［7］Huang L，So H C.Source enumeration via MDL criterion based on linear shrinkage estimation of noise subspace covariance matrix［J］.IEEE Transactions on Signal Processing，2013，61（19）：4806－4821.

［8］Lu Z H，Zoubir A M.Source enumeration using the PDF of sample eigenvalues via information theoretic criteria［C］∥IEEE International Conference on Acoustics Speech Signal Processing，Kyoto，Japan，2012，2：3361－3364.

［9］Pesavento M，Gershman A B.Maximum－likelihood direction－of－arrival estimation in the presence of unknown nonuniform noise［J］.IEEE Transactions on Signal Processing，2001，49（7）：1310－1324.

［10］Wu H T，Yung J F，Chen F K.Source number estimator using gerschgorin disks［J］.International Conference on Acoustics Speech Signal Processing，Adelaide Australia，1994，4（4）：261－264.

［11］Wu H T，Yung J F，Chen F K.Source number estimator using transformed gerschgorin radii［J］.IEEE Transactions On Signal Processing，1995，43（6）：1325－1333.

［12］Yung J F，Wu H T.Gerschgorin radii based source number detection for closely spaced signals［C］∥IEEE Transactions Acoustics Speech Signal Processing，Atlanta May，1996：3054－3057.

［13］Gu J F，Wei P，Tai H M.Detection of the number of sources at low signal－to－noise ratio［C］∥IET Signal Processing，2007，1（1）：2－8.

［14］Zoubir A，Boashash B.The bootstrap and its application in signal processing［C］∥IEEE Signal Processing Magazine，1998，1：55－76.

［15］Politis D.Computer－intensive methods in statistical analysis［J］.IEEE Signal Processing Magazine，1998，1：39－55.

［16］Aouada S，Zoubir A M，See C M S.Source detection in the presence of nonuniform noise［C］∥IEEE International Conference on Acoustics Speech Signal Processing，Montreal，Canada，2004，2：165－168.

［17］Scharf L L.Statistical Signal Processing：Detection，Estimation and Time Series Analysis［M］.Addison Wesley，1990.

［18］Zoubir A，Iskander D.Bootstrap Techniques for Signal Processing［M］.Cambridge：Cambridge University Press，2004.

［19］Aouada S，Transkov D，Heureuse N，et al.Application of the bootstrap to source detection in Nonuniform noise［C］∥IEEE International Conference on Acoustics Speech Signal Processing，Pennsylvania，USA，2005，4：997－1000.

［20］Brcich R，Zoubir A M，Pelin P.Detection of sources using bootstrap techqiques［J］.IEEE Transactions on Signal Processing，2002，50（2）：206－215.

［21］Eforn B，Tibshirani R.An Introduction to the Bootstrap［M］.London：Chapman＆Hall，1993.