夏文杰 蔡志明

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

1 引言

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

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

2 功率譜熵

2.1 定義

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

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

(1)

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

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

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

整個(gè)檢測(cè)過(guò)程可看為一個(gè)二元假設(shè)檢驗(yàn)問(wèn)題：

H0:x(n)=w(n)

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

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

(2)

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

針對(duì)正弦信號(hào)，s(n)可表示為：

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

(3)

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

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

(4)

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

(5)

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

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

(6)

其中，

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

(7)

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

(8)

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

3 功率譜熵檢測(cè)的概率密度分布

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

3.1 零假設(shè)

H0假設(shè)下有：

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

(9)

k=0，1，…，N-1

(10)

式中，

(11)

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

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

(12)

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

(13)

(14)

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

(15)

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

對(duì)于Yk概率密度函數(shù)可用隨機(jī)變量的函數(shù)中的定義法求出。利用定義法求Yk的概率密度函數(shù)時(shí)，需要用到y(tǒng)=xlnx的反函數(shù)。但y=xlnx在其定義域內(nèi)不具有單調(diào)性：在x∈(0，1/e)時(shí)，y單調(diào)遞減；在x∈(1/e，+)時(shí)，y單調(diào)遞增。y的最小值在x=e-1時(shí)取到，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)時(shí)，y=xlnx的反函數(shù)為：

(16)

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

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

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

(17)

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

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

(18)

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

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

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

(19)

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

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

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

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

(20)

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

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

(21)

(22)

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

μ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)性的原因有兩個(gè)方面：

(24)

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

(25)

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

(26)

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

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

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

圖2 檢驗(yàn)統(tǒng)計(jì)量T均值和方差理論計(jì)算結(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)為信號(hào)的能量全部集中在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)計(jì)的方法得到。Yk′和Yk″的均值、方差可根據(jù)相應(yīng)的概率密度函數(shù)求出。觀察式(29)、(31)可知，當(dāng)信噪比較低時(shí)，兩式形式相差不大，檢測(cè)統(tǒng)計(jì)量T在H1假設(shè)下仍可近似認(rèn)為服從高斯分布，顯然，當(dāng)信噪比較高時(shí)，k″在頻譜中占主要成分，不滿足中心極限定理，T的分布需進(jìn)一步分析。

4 功率譜熵檢測(cè)接收機(jī)工作特性曲線

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

(33)

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

(34)

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

(35)

檢測(cè)概率為

(36)

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

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

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

圖3 功率譜熵檢測(cè)性能分析

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

圖4 功率譜熵檢測(cè)方法和能量檢測(cè)方法性能對(duì)比圖

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

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

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

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

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

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

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

5 結(jié)論

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

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