夏文杰 蔡志明

(海軍工程大學(xué)水聲工程教研室，湖北武漢 430033)

1 引言

功率譜熵作為表征信號復(fù)雜度的非線性特征量，已經(jīng)用于運動物體的檢測[1]、轉(zhuǎn)子振動故障分析[2]、無線電頻譜感知問題[3]、語音信號端點檢測[4]、診斷氣液兩相流流型變化[5]、水聲脈沖信號檢測[6]、船舶軸頻電場線譜提取[7]等信號處理領(lǐng)域。對于噪聲而言，混亂程度高，其功率譜熵大；對于信號而言，混亂程度低，其功率譜熵小。該方法利用信號與噪聲之間功率譜熵大小的差異，判斷是否存在信號。通過對接收數(shù)據(jù)進行離散傅里葉變換得到功率譜，對其歸一化后再做功率譜熵分析，并以此作為檢測統(tǒng)計量。此方法可在無先驗信息及低信噪比條件下，實現(xiàn)單頻或調(diào)頻信號脈沖的檢測。

本文針對白高斯噪聲中的正弦信號，從理論上對功率譜熵檢測方法進行分析，計算得到零假設(shè)H0和備擇假設(shè)H1下檢測統(tǒng)計量的概率密度函數(shù)，并給出虛警概率Pfa，檢測門限γ、檢測概率Pd的相應(yīng)表達式。根據(jù)紐曼皮爾遜準(zhǔn)則，在給定虛警概率條件下，可求出相應(yīng)的檢測門限以及檢測概率。將檢測門限和檢測概率的理論計算結(jié)果和蒙特卡洛實驗仿真結(jié)果進行對比，吻合度較好，驗證了理論的合理性和準(zhǔn)確性。當(dāng)Pfa=10-4、SNR=-13 dB時，Pd大約等于0.5，該方法檢測性能要優(yōu)于能量檢測。對海試數(shù)據(jù)進行處理分析，通過該方法能在低信噪比條件下有效檢測信號。

2 功率譜熵

2.1 定義

Shannon把熵的概念引入信息論中，表示信息系統(tǒng)中描述信息的能力。一個系統(tǒng)有序程度越高，熵就越小，所含的信息量就越大；反之，系統(tǒng)的無序?qū)?yīng)熵值大。信息對應(yīng)不確定性，而不確定性在概率論中是用隨機事件或隨機變量來描述的。

X是取值為有限的隨機變量，其概率分布律為P{X=xi}=pi，i=1，2，…，n。熵定義為：

(1)

式中，a為對數(shù)底數(shù)，可取2、e、10。本文取為e以便數(shù)學(xué)演繹。

使用Shannon信息熵概念定量地計算信號功率譜的不確定性，即從頻域角度定義的信息熵，可度量被測信號在頻域上的復(fù)雜程度，這就是所謂的功率譜熵。利用功率譜熵表征和提取信號與背景噪聲的不確定性差異。具體地，若被測數(shù)據(jù)中不含信號，其功率譜熵值就應(yīng)較大；否則功率譜熵值較小。在實際中，被測數(shù)據(jù)總是有限長的，必然存在柵欄效應(yīng)從而導(dǎo)致譜泄露。當(dāng)譜泄露時，該方法的性能有所下降，但仍然在可接受的范圍內(nèi)。為簡便這里假設(shè)信號未發(fā)生譜泄露。

2.2 檢測統(tǒng)計量構(gòu)造

整個檢測過程可看為一個二元假設(shè)檢驗問題：

H0:x(n)=w(n)

H1:x(n)=s(n)+w(n)

n=0，1，2，…，N-1

(2)

式中，x(n)表示觀測值，s(n)表示信號，w(n)為加性噪聲，N為樣本長度，取偶數(shù)。

針對正弦信號，s(n)可表示為：

s(n)=Acos(ω0n+φ)

(3)

式中所有參數(shù)均未知，A為幅度，ω0=2πf0/fs為歸一化角頻率, f0和fs分別為信號頻率和采樣率，φ為相位?；诠β首V熵檢測的步驟如下：

(1)利用離散傅里葉變換得到信號功率譜密度的估計：

(4)

式中，X(k)為信號的離散傅里葉變換。

(5)

式中，p(k)表示第k個功率譜在總功率譜中占的比值。

(3)計算出相應(yīng)的功率譜信息熵，簡稱功率譜熵H：

(6)

其中，

Yk=p(k)ln[p(k)]

(7)

(4)H作為檢測統(tǒng)計量T：

(8)

式中，γ為檢測門限。當(dāng)T[x(n)]大于γ時，則認(rèn)為接收數(shù)據(jù)中不存在正弦信號；當(dāng)T[x(n)]小于γ時，則認(rèn)為接收數(shù)據(jù)中存在正弦信號。

3 功率譜熵檢測的概率密度分布

根據(jù)功率譜熵的定義，檢測統(tǒng)計量T[x(n)]是由每個頻點的功率譜值經(jīng)過若干變換得到Y(jié)k后再相加而來。Yk是一個隨機變量，其兩兩之間具有一定的相關(guān)性，但是隨著采樣點數(shù)N的增加，兩兩之間的相關(guān)性將減弱。當(dāng)N足夠大(實際上只要N達到256)時，Yk之間足以解除因歸一化產(chǎn)生的線性相關(guān)，根據(jù)中心極限定理，檢測統(tǒng)計量T[x(n)]將近似服從正態(tài)分布。零假設(shè)和備擇假設(shè)的數(shù)字特征有所不同，且由于相關(guān)性的原因?qū)τ趦煞N假設(shè)的方差需要進行修正，將在下面進行討論。

3.1 零假設(shè)

H0假設(shè)下有：

x(n)=w(n)，n=0，1，2，…，N-1

(9)

k=0，1，…，N-1

(10)

式中，

(11)

對于白噪聲序列，每個點之間都是相互獨立的隨機變量，因此w(n)的離散傅里葉變換W(k)對于每一個離散頻率k而言，可以看作是相互獨立的隨機變量的線性組合，其結(jié)果仍為隨機變量[8-9]。因為高斯分布的線性組合仍然是高斯分布，所以W(k)，WR(k)，WI(k)也服從高斯分布。WR(k)和WI(k)的均值、方差為[8]：

E[WR(k)]=E[WI(k)]=0

(12)

w(n)的功率譜估計為：

(13)

(14)

根據(jù)檢驗統(tǒng)計量的構(gòu)造，還須將PX(k)進行歸一化：

(15)

可知，p(k)服從參數(shù)為N的指數(shù)分布，即p(k)～Exp(N)。

對于Yk概率密度函數(shù)可用隨機變量的函數(shù)中的定義法求出。利用定義法求Yk的概率密度函數(shù)時，需要用到y(tǒng)=xlnx的反函數(shù)。但y=xlnx在其定義域內(nèi)不具有單調(diào)性：在x∈(0，1/e)時，y單調(diào)遞減；在x∈(1/e，+)時，y單調(diào)遞增。y的最小值在x=e-1時取到，ymin=-e-1。根據(jù)反函數(shù)的定義y=xlnx(x>0)不存在反函數(shù)?？紤]到y(tǒng)=xlnx是分段單調(diào)的，這里把y=xlnx分為兩段，再分別在其單調(diào)區(qū)間內(nèi)求反函數(shù)。

在區(qū)間(0，1/e)時，y=xlnx的反函數(shù)為：

(16)

其中x+lnx=y的解為x=wrightOmega(y)，為便于書寫，記x=s1(y)。

當(dāng)y<0時，根據(jù)復(fù)對數(shù)運算法則：

lny=ln|y|+j(arg(y)+2πk)，k∈Z

(17)

式中，Z表示整數(shù)。這將導(dǎo)致多值性的出現(xiàn)，不滿足函數(shù)唯一性的要求。解決復(fù)對數(shù)多值性的一般辦法是取主值進行運算，將arg(y)對π求模得到主值，用ARG(y)表示，ARG(y)∈(-π，π)。上式化為：

lny=ln|y|+jARG(y)，k∈Z

(18)

此時滿足函數(shù)值唯一性的性質(zhì)。

同理可求在區(qū)間(1/e，+)時，y=xlnx的反函數(shù)為：

x=exp(lambertw(0， y))， y∈(-1/e，+)

(19)

其中xex=y的解為x=lambertw(0， y)，記x=s2(y)。

根據(jù)隨機變量的函數(shù)分布中的定義法可得求Yk的分布函數(shù)和概率密度函數(shù)：

(1)Yk的分布函數(shù)

由定義法和數(shù)形結(jié)合法可得Y的分布函數(shù)：

(20)

(2)Yk的概率密度函數(shù)

由變限積分的求導(dǎo)公式對分布函數(shù)FYk(y)求導(dǎo)可得Yk的概率密度函數(shù)：

(21)

(22)

根據(jù)隨機變量Yk的概率密度函數(shù)，分段積分可求出其均值E[Y0]和方差D[Y0]。檢驗統(tǒng)計量的均值和方差為：

μ0=E[T]=-NE[Y0]

D[2(Y0+Y1+…+YN/2-1)]=

2ND[Y0](1+(N/2-1)ρYiYj)

(23)

式中， ρYiYj為Yi、Yj間的相關(guān)系數(shù)(i≠j)。Yk之間產(chǎn)生相關(guān)性的原因有兩個方面：

(24)

式中，cov(Yi，Yj)為Yi與Yj的協(xié)方差(i+j=N-1)。為消除此相關(guān)性，取p(k)進行功率譜統(tǒng)計時，僅取前一半。

(25)

這里，M為除p(i)外其他所有隨機變量(j≠i;i, j≤N/2-1)之和，可得

(26)

圖1 隨機變量之間的相關(guān)系數(shù)

表1給出了一些常用采樣點數(shù)的Yk之間的相關(guān)系數(shù)。

表1 Yk之間的相關(guān)系數(shù)

圖2 檢驗統(tǒng)計量T均值和方差理論計算結(jié)果和蒙特卡洛仿真結(jié)果比較

3.2 備擇假設(shè)

H1假設(shè)下有：

x(n)=s(n)+w(n)，n=0，1，2，…，N-1

(27)

=

WR(k)-jWI(k)，k=0，1，…，N-1

(28)

為了便于分析，假設(shè)沒有譜泄露，認(rèn)為信號的能量全部集中在k=Nω0/2π處。

fp(k′)(y)=N(1+SNR)exp(-N(1+SNR)y)，y>0

(29)

(30)

(31)

μ1=E[T]=-(N-2)E[Yk′]-2E[Yk″]

D[2(Y0+Y1+…+YN/2-1)]=

4[(N/2-1)D[Yk′]+D[Yk″]+

(N-2)(2+(N-4)ρk′)D[Yk″]+

(32)

式中， ρk′為Yk′間的相關(guān)系數(shù)，ρk″為Yk′與Yk″的相關(guān)系數(shù)。與3.1節(jié)中一樣，相關(guān)系數(shù)用統(tǒng)計的方法得到。Yk′和Yk″的均值、方差可根據(jù)相應(yīng)的概率密度函數(shù)求出。觀察式(29)、(31)可知，當(dāng)信噪比較低時，兩式形式相差不大，檢測統(tǒng)計量T在H1假設(shè)下仍可近似認(rèn)為服從高斯分布，顯然，當(dāng)信噪比較高時，k″在頻譜中占主要成分，不滿足中心極限定理，T的分布需進一步分析。

4 功率譜熵檢測接收機工作特性曲線

根據(jù)3.1，3.2節(jié)的分析，檢測統(tǒng)計量T服從如下分布：

(33)

采用紐曼皮爾遜(N-P)準(zhǔn)則實現(xiàn)假設(shè)檢驗。首先確定虛警概率Pfa，然后計算不同信噪比情況下的檢測概率Pd來衡量檢測性能。通過給定的虛警概率Pfa可以確定出檢測門限γ。根據(jù)Pfa=P(T<γ|H0)，得

(34)

γ=σ0Q-1(Pfa)+μ0

(35)

檢測概率為

(36)

將式(34)得到的檢測門限γ代入式(35)中可以得檢測概率Pd。

4.1 仿真數(shù)據(jù)驗證

仿真實驗具體參數(shù)如下：信號為正弦信號，噪聲為高斯白噪聲，采樣點數(shù)N=1024，采樣頻率fs=1024Hz。圖3(a)給出了虛警概率Pfa從10-4增加到10-3時，理論計算得到的檢測門限γ；圖3(b)給出了在不同信噪比條件下功率譜熵檢測性能，信噪比由-20dB以-2dB的間隔增加至-14dB，理論計算結(jié)果和蒙特卡洛仿真結(jié)果吻合較好，驗證了模型的準(zhǔn)確性；圖3(c)給出了不同虛警概率條件下功率譜熵檢測性能。從圖中可以看出，當(dāng)信噪比超過-10dB時，即便十萬分之一的虛警水平，其檢測概率可接近于1。

圖3 功率譜熵檢測性能分析

圖4給出了能量檢測[11]與功率譜熵檢測性能。仿真的各參數(shù)與之前一致，設(shè)置虛警概率Pfa=10-4。對比可知，檢測概率Pd=0.5時要求功率譜熵檢測方法的信噪比為-13dB，比能量檢測方法優(yōu)4dB。

圖4 功率譜熵檢測方法和能量檢測方法性能對比圖

4.2 海試數(shù)據(jù)檢驗

海試實驗：聲源船發(fā)射頻率為500Hz的CW脈沖信號，信號出現(xiàn)的時間為42s至48s，脈寬6s。數(shù)據(jù)時長為1min，采樣頻率fs=8000Hz，工作頻帶為400Hz～800Hz。

圖5是常規(guī)預(yù)成波束的輸出，在131°方位明顯存在一段的信號，其出現(xiàn)的時間大約為42s至48s之間。用此方位的波束數(shù)據(jù)進行功率譜熵檢測的自動處理。

圖6是在131°方位波束域數(shù)據(jù)做功率譜熵處理的結(jié)果。在400Hz～800Hz頻段內(nèi)，背景噪聲不符合白噪聲的條件，理論模型對此情形有所失配。為此，選取一段無CW脈沖也無強干擾的背景數(shù)據(jù)，通過數(shù)據(jù)學(xué)習(xí)確定在虛警概率Pfa=10-4時的檢測門限γ=5.013(理論檢測門限γ=5.427)。從圖中可認(rèn)為在時間42.01s至48.57s時存在信號，與現(xiàn)實情況相符。若直接以此來估計脈寬，可粗略地判斷脈寬為6.56s，精度為9.3%。

圖5 遠(yuǎn)場常規(guī)波束形成

圖6 功率譜熵處理結(jié)果

圖7 海試數(shù)據(jù)和仿真數(shù)據(jù)檢測性能對比

5 結(jié)論

本文基于功率譜熵檢測方法，針對白高斯噪聲背景中檢測未知正弦信號的問題，從理論上推導(dǎo)了零假設(shè)和備擇假設(shè)下檢測統(tǒng)計量的概率密度函數(shù)，給出虛警概率Pfa、檢測門限γ和檢測概率Pd相應(yīng)的表達式。與蒙特卡洛實驗仿真結(jié)果比較，吻合度高，驗證了理論結(jié)果的準(zhǔn)確性。對信噪比較高的情況下，備擇假設(shè)所求得的概率密度函數(shù)不夠準(zhǔn)確，需要進行進一步分析處理。將功率譜熵檢測方法與能量檢測方法的性能進行比較，通過仿真驗證可以看出該方法在性能上要優(yōu)于后者，性能相差4 dB。

對海上試驗數(shù)據(jù)進行處理分析。通過對背景噪聲學(xué)習(xí)獲得的檢測門限與理論門限接近，其偏差主要源于背景有色功率譜相對白譜的差別。用數(shù)據(jù)自學(xué)習(xí)的門限值進行自動檢測處理，在虛警概率Pfa=10-4，檢測概率Pd=0.5時，要求信噪比為-10 dB。驗證了功率譜熵檢測方法在實際應(yīng)用場景中的可檢測信噪比與理論預(yù)報不超過2 dB。