α穩(wěn)定分布特征指數(shù)估計(jì)算法

關(guān)濟(jì)實(shí),石要武,邱建文,單澤彪,史紅偉

(1.吉林大學(xué) 通信工程學(xué)院,長春130022;2.中廣核研究院有限公司 北京分公司 北京100083;3.長春理工大學(xué) 電子信息工程學(xué)院,130022)

0 引 言

一類帶有沖擊特性的非高斯噪聲在大氣、水下聲學(xué)以及電磁干擾等噪聲中大量存在,該類噪聲表現(xiàn)出極強(qiáng)的脈沖特性,由于α穩(wěn)定分布可以用來很好地描述這一類噪聲[1-3],近年來α穩(wěn)定分布噪聲背景下的信號(hào)處理方法得到了廣泛研究。在α噪聲背景下的信號(hào)處理方法的研究中,分?jǐn)?shù)低階矩和共變被大量應(yīng)用[4-6],而分?jǐn)?shù)低階矩和共變只有在α穩(wěn)定分布噪聲參數(shù)已知的情況下才有意義,所以大多研究都假設(shè)α噪聲的參數(shù)已知,至少是其特征指數(shù)α已知。在實(shí)際應(yīng)用中,由于沒有α噪聲的相關(guān)參數(shù),導(dǎo)致這些算法的應(yīng)用受到很大的限制,因此研究α穩(wěn)定分布的參數(shù)估計(jì)方法對于以上這些算法的應(yīng)用具有支撐作用。在α穩(wěn)定分布參數(shù)估計(jì)方面,Tsihrintzis[7]、Ma等[8]基于分?jǐn)?shù)低階矩和負(fù)階矩理論提出了多種通過觀測序列估計(jì)α穩(wěn)定參數(shù)的方法,Kuruoglu[9]借助混合高斯模型,提出了α穩(wěn)定概率密度函數(shù)的一種新的表達(dá)式,并利用該表達(dá)式進(jìn)行α穩(wěn)定分布的參數(shù)估計(jì)。Du Mouchel[10]提出了一種最大似然估計(jì)方法,Brorsen等[11]對該方法進(jìn)行蒙特卡羅仿真,獲得了相當(dāng)好的結(jié)果,然而,該方法屬于高度的非線性優(yōu)化問題,由于需要計(jì)算復(fù)雜的數(shù)值積分,計(jì)算量非常大。本文提出一種α穩(wěn)定分布信號(hào)的α和γ參數(shù)估計(jì)的方法,該方法利用α的疊加性質(zhì),通過簡單的運(yùn)算即可實(shí)現(xiàn)對α估計(jì),容易理解且計(jì)算量小,估計(jì)精度高。

1 α穩(wěn)定分布的定義和性質(zhì)

1.1 α穩(wěn)定分布的定義

α穩(wěn)定分布的概念是Levy在1925年研究廣義中心極限定理時(shí)提出的[3],中心極限定理引出了高斯分布,而廣義中心極限定理則是放開了中心極限定理中有限方差這一限制后形成的。α穩(wěn)定分布是負(fù)責(zé)廣義中心極限定理的一簇分布,高斯分布是α穩(wěn)定分布在α＝2時(shí)的一個(gè)特例[12]。除了有限的幾種情況,α穩(wěn)定分布沒有封閉的概率密度函數(shù)表達(dá)式[13],α穩(wěn)定分布的定義是由其概率密度函數(shù)的傅里葉變換,即特征函數(shù)給出的。

定義 若一個(gè)隨機(jī)變量X的特征函數(shù)可以表示為:

式中:－∞＜μ＜∞,γ＞0,0＜α≤2,－1≤β≤1,這4個(gè)參數(shù)決定了該分布的形狀,sign(x)為符號(hào)函數(shù)。稱該變量X服從α穩(wěn)定分布,服從α穩(wěn)定分布的噪聲稱為沖擊噪聲,脈沖噪聲或穩(wěn)定分布噪聲,寫作:

式中:α稱為特征指數(shù),它決定了脈沖特性的程度,α越小,脈沖程度越高;β稱為對稱系數(shù),當(dāng)β＜0時(shí),信號(hào)的概率密度函數(shù)會(huì)出現(xiàn)左傾,反之出現(xiàn)右傾。當(dāng)β＝0時(shí),α穩(wěn)定分布的概率密度函數(shù)呈現(xiàn)左右對稱形狀,稱為對稱α穩(wěn)定分布(SαS);μ稱為位置參數(shù),表示α穩(wěn)定分布概率密度函數(shù)的中心,也是α穩(wěn)定分布變量的期望值;γ稱為分散系數(shù),在同一個(gè)α值下,γ越大,則概率密度函數(shù)的兩側(cè)延伸越嚴(yán)重,中間部分的值則會(huì)相對降低,與高斯噪聲的方差類似。

1.2 α穩(wěn)定分布的性質(zhì)

現(xiàn)有文獻(xiàn)對α穩(wěn)定分布的性質(zhì)做了詳細(xì)的研究,本文不再一一說明,僅介紹本文算法相關(guān)的兩條性質(zhì),即疊加性質(zhì)和p范數(shù)性質(zhì)[14]。

性質(zhì)1 疊加性質(zhì),若X1～Sα(μ1,β1,γ1),X2～Sα(μ2,β2,γ2),則它們的和變量X＝X1＋X2服從X～Sα(μ,β,γ),其中:

性質(zhì)2 若X～Sα(u,β,γ),0＜α＜2且當(dāng)α＝1時(shí),有β＝0,則對于所有0＜p＜α,有:

2 基于組合分布和分?jǐn)?shù)低階矩的系數(shù)估計(jì)方法

2.1 特征指數(shù)估計(jì)算法原理

根據(jù)α分布的性質(zhì)1,兩個(gè)獨(dú)立的同特征指數(shù)的穩(wěn)定分布變量相加時(shí),其和變量同樣服從α分布,并且其分散系數(shù)仍為α,當(dāng)多個(gè)同分布的α信號(hào)相加時(shí),這一性質(zhì)可以以如下形式表述:

推理 當(dāng)X1,X2,…,X N相互獨(dú)立且均服從于同一分布Sα(μ,β,γ)分布時(shí),其和分布X＝X1＋X2＋…＋X N服從Sα(μ0,β0,γ0),則其參數(shù)表達(dá)式為:

對照性質(zhì)1,此推理很容易得證。

從以上推理中可以看到,多個(gè)獨(dú)立同分布的α變量的和變量仍服從α分布,且其對稱系數(shù)不變。注意到這一性質(zhì)分散系數(shù)的表達(dá)式中并沒有對稱系數(shù),所以分散系數(shù)的這一性質(zhì)對于非對稱α分布也是適用的。

綜合以上兩性質(zhì)和推理,多個(gè)同分布變量相加時(shí)有:

式中:i＝1,2,…,N。

和變量X與原變量X i的p階矩之間具有以下關(guān)系:

由式(13)可知,當(dāng)有多個(gè)同分布α穩(wěn)定分布隨機(jī)變量時(shí),可以將多個(gè)隨機(jī)變量組合,利用組合變量與原變量之間p范數(shù)關(guān)系,得到特征指數(shù)α的估計(jì)。

2.2 α穩(wěn)定分布參數(shù)估計(jì)方法

2.2.1 參數(shù)估計(jì)方法

式(13)表明多個(gè)同分布的α信號(hào)相加時(shí),其統(tǒng)計(jì)量E[|X|p]1/p之間的比只與參與相加的變量個(gè)數(shù)和特征指數(shù)α有關(guān)。利用這一特征,可以從實(shí)際采樣信號(hào)中估計(jì)得到特征指數(shù)α。

當(dāng)一個(gè)變量的時(shí)間相關(guān)性很低時(shí),相隔足夠長時(shí)間以外的多個(gè)采樣值之間可以認(rèn)為是不相關(guān)的。對于一個(gè)α分布的變量,相隔足夠長的時(shí)間進(jìn)行多次采樣,則可以認(rèn)為多次采樣結(jié)果為多個(gè)獨(dú)立同分布的α分布變量,這多次采樣相加,得到的新變量與原變量之間可以根據(jù)式(13)得到α的估計(jì)。對于一個(gè)足夠長時(shí)間的采樣序列X,該方法可以描述如下:

設(shè)有序列長度為L的α穩(wěn)定分布變量采樣序列X＝{x i,i＝1,2,…,L},將X等分為N＋1份,每份長度為K(K＝L/(N＋1)),則有:

由以上定義可知,Y和Z是相互獨(dú)立的隨機(jī)變量。

算法的執(zhí)行步驟如下:

(1)對待測信號(hào)采樣L次,形成X。

(2)將L次信號(hào)平均分割為N段,形成X m,m＝1,2,…,N。

(3)將分割的數(shù)據(jù)相加,形成X′。

(4)對X和X′分別計(jì)算統(tǒng)計(jì)量T和T′。

(5)將數(shù)據(jù)代入式(19)得到特征指數(shù)α的估計(jì)值。

該方法不需要復(fù)雜的計(jì)算,相比于現(xiàn)有的參數(shù)估計(jì)方法,計(jì)算量較小。并且不涉及α分布概率密度函數(shù)的封閉形式的近似,原理簡單可靠,便于理解和編程實(shí)現(xiàn)。在實(shí)際執(zhí)行時(shí),由于α未知,為了使p階矩收斂,可以使p盡量小。

2.2.2 估計(jì)方法的無偏性與一致性分析

式(12)由α穩(wěn)定分布的性質(zhì)導(dǎo)出,因此它是嚴(yán)格成立的,但在實(shí)際運(yùn)算中,無法直接得到變量的p階矩,只能以式(17)得到的估計(jì),式(17)中統(tǒng)計(jì)量T′的數(shù)學(xué)期望為:

再來分析估計(jì)的一致性,定義變量D如下:

α與D之間是冪函數(shù)關(guān)系,二者之間是一一對應(yīng)的,因此如果算法中對D的估計(jì)是無偏的,則對α的估計(jì)也將是無偏的,算法中D的估計(jì)值為:

當(dāng)2p＜α?xí)r,均為有限值,當(dāng)采樣長度無窮長,即L無窮大時(shí),K＝L/N也趨于無窮大,因此式(28)可寫作:

因此,當(dāng)采樣長度無窮大時(shí),式(27)可表示為:

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

實(shí)驗(yàn)1 不同α值的特征指數(shù)估計(jì)實(shí)驗(yàn)

實(shí)驗(yàn)?zāi)康?驗(yàn)證算法的有效性。

實(shí)驗(yàn)方法:α的取值從0.1至2,每隔0.1取一個(gè)α值,利用數(shù)值方法實(shí)現(xiàn)β＝0,γ＝1,長度為10 000的α分布序列,分割為兩段,通過式(19),取p＝0.1計(jì)算α的估計(jì)值,并與產(chǎn)生數(shù)據(jù)所用的α相比較,以驗(yàn)證算法的有效性。

估計(jì)結(jié)果與生成數(shù)據(jù)時(shí)所用α之間關(guān)系如圖1所示。

從圖1可以看到,該方法準(zhǔn)確地估計(jì)了α的值。為考察算法精度,定義估計(jì)均方誤差:

圖1 基本的α估計(jì)實(shí)驗(yàn)結(jié)果Fig.1 Estimation result ofαwith basic experiment

每個(gè)α取值進(jìn)行100次蒙特卡羅實(shí)驗(yàn),計(jì)算均方誤差,圖2為均方誤差圖。

圖2 α估計(jì)均方誤差Fig.2 Estimation MSE ofαwith basic experiment

可以看到,隨著α變大,均方誤差逐漸變大,當(dāng)α＝2時(shí),均方誤差僅為－12 dB,具有較高的估計(jì)精度。

實(shí)驗(yàn)2 分割段數(shù)對算法的影響

實(shí)驗(yàn)?zāi)康?驗(yàn)證分割段數(shù)對于算法精度的影響。

實(shí)驗(yàn)方法:特征指數(shù)取0.6,0.9,1.2,1.6,分段數(shù)分別取[2,3,4,5,6,8,10,12,15,16,20,24]進(jìn)行計(jì)算,每種情況進(jìn)行100次蒙特卡羅實(shí)驗(yàn),為了得到完整的分段,實(shí)驗(yàn)數(shù)據(jù)長度取12 000,記錄估計(jì)值和均方誤差,觀察分段數(shù)對于估計(jì)誤差的影響。

圖3為分段數(shù)變化時(shí)特征指數(shù)的估計(jì)結(jié)果和均方誤差。

從圖3可以看出,不同分段數(shù)時(shí),本方法均可以有效估計(jì)α值。從均方誤差結(jié)果圖中可以看到,段數(shù)變大估計(jì)均方誤差變大,這是由于分段數(shù)變大使得每段數(shù)據(jù)量變小,從而使得p范數(shù)的估計(jì)誤差變大引起的。本方法對α的估計(jì)誤差始終在－13 dB以下,估計(jì)精度較高。

實(shí)驗(yàn)3 分散系數(shù)對算法的影響實(shí)驗(yàn)

實(shí)驗(yàn)?zāi)康?驗(yàn)證不同分散系數(shù)γ時(shí)的算法性能。

實(shí)驗(yàn)方法:γ從0.1變到10,每隔0.1取一個(gè)點(diǎn),α取0.6,0.9,1.2,1.6,進(jìn)行4次實(shí)驗(yàn),每種情況進(jìn)行100次蒙特卡羅實(shí)驗(yàn),分段數(shù)取4,每次實(shí)驗(yàn)數(shù)據(jù)長度為10 000,記錄估計(jì)值和均方誤差并形成曲線。

圖4為γ分段數(shù)變化時(shí)估計(jì)誤差結(jié)果。

圖3 分段數(shù)對α估計(jì)結(jié)果的影響Fig.3 Influence of section number on estimation ofα

圖4 分散系數(shù)對α估計(jì)結(jié)果的影響Fig.4 Influence ofγon estimation ofα

從圖4可以看到,γ從0.1到10,本方法始終能正確估計(jì)α值。從估計(jì)均方誤差結(jié)果可以看出,不同γ值下,α值估計(jì)的均方誤差變化不大,且始終能保持－14 d B以下的均方誤差,估計(jì)精度較高。

實(shí)驗(yàn)4 高斯噪聲與SαS噪聲共存時(shí)的估計(jì)效果

實(shí)驗(yàn)?zāi)康?驗(yàn)證高斯信號(hào)與SαS信號(hào)共存時(shí)算法對于α的估計(jì)效果。

實(shí)驗(yàn)方法:α取0.6,0.9,1.2,1.6,在SαS信號(hào)中摻雜入廣義信噪比不同的高斯噪聲,其他條件與前述實(shí)驗(yàn)相同,依據(jù)本算法估計(jì)信號(hào)α值。每種情況進(jìn)行100次蒙特卡羅實(shí)驗(yàn),記錄α的估計(jì)結(jié)果均值和均方誤差。

圖5為高斯噪聲摻雜時(shí)α的估計(jì)結(jié)果。

圖5 高斯噪聲摻雜時(shí)算法對α的估計(jì)結(jié)果Fig.5 Estimation ofαwhen Gaussian noise added

從圖5可以看出,在信噪比0 d B以下,高斯噪聲的添加對特征指數(shù)α估計(jì)結(jié)果的影響不大,估計(jì)均方誤差在－8 d B以下。說明本算法可以估計(jì)高斯隨機(jī)變量與α穩(wěn)定分布隨機(jī)變量共存時(shí)α穩(wěn)定分布隨機(jī)變量的特征指數(shù)。

4 結(jié)束語

針對α穩(wěn)定分布信號(hào)處理中的參數(shù)估計(jì)問題,本文提出了基于α穩(wěn)定分布p-范數(shù)和疊加性質(zhì)的特征指數(shù)α估計(jì)方法。將α穩(wěn)定分布隨機(jī)變量的采集序列等分相加,結(jié)果序列的p范數(shù)與原序列的p范數(shù)之間存在的固定關(guān)系可以用來估計(jì)特征指數(shù)α。該方法不需要β的先驗(yàn)知識(shí)即可估計(jì)α值,在很大范圍內(nèi)可以正確估計(jì)α值。在0.1～2范圍內(nèi),α值估計(jì)的均方誤差均小于－13 d B,分段數(shù)和隨機(jī)變量的γ變化對估計(jì)性能影響不大,在一定程度的高斯噪聲與α穩(wěn)定分布摻雜時(shí),仍能正確估計(jì)α值。本文詳細(xì)說明了算法原理和算法步驟,并證明了估計(jì)算法的無偏性和一致性。該方法可以為針對α噪聲的信號(hào)處理算法提供特征指數(shù)α的先驗(yàn)知識(shí),為基于分?jǐn)?shù)低階矩的算法提供階次選擇的依據(jù)。

[1]Schilder M.Some structure theorems for the symmetric stable laws[J].The Annals of Mathematical Statistics,1970,41(2):412-421.

[2]Weron A.Stable Processes and Measures:a Survey[M].Berlin:Springer,1983:306-364.

[3]Shao M,Nikias C L.Signal processing with fractional lower order moments:stable processes and their applications[J].Proceeding of the IEEE,1993,81(7):986-1010.

[4]孫永梅.α穩(wěn)定分布參數(shù)估計(jì)與譜分析理論及應(yīng)用研究[D].大連:大連理工大學(xué)信息與通信工程學(xué)院,2006.Sun Yong-mei.Study on the theory and application of parameter estimation and spectral analysis ofα-stable distribution[D].Dalian:School of Information and Communication Engineering,Dalian University of Technology,2006.

[5]李森.穩(wěn)定分布噪聲下通信信號(hào)處理新算法及性能分析[D].大連:大連理工大學(xué)信息與通信工程學(xué)院,2011.Li Sen.Communication signal processing new algorithms and performance analysis in stable noise[D].Dalian:School of Information and Communication Engineering,Dalian University of Technology,2011.

[6]何勁.α穩(wěn)定分布噪聲背景下陣列信號(hào)處理方法研究[D].南京:南京理工大學(xué)電子工程與光電技術(shù)學(xué)院,2007.He Jin.Investigation on array signal processing a midα-stable noise[D].Nanjing:School of Electronic and Optical Engineering,Nanjing University of Science and Technology,2007.

[7]Tsihrintzis G A,Nikias C L.Fast estimation of the parameters of alpha-stable impulsive interference[J].IEEE Trans on Signal Processing,1996,44(6):1492-1503.

[8]Ma X,Nikias C L.Parameter estimation and blind channel identification in impulsive interference[J].IEEE Trans on Signal Processing,1995,43(11):2884-2897.

[9]Kuruoglu E E.Analytical representation for positiveα-stable densities[C]∥IEEE International Conference on Acoustics,Hong Kong,2003:VI-729-32.

[10]Du Mouchel W H.Stable distributions in statistical inference:information from stably distributed samples[J].Journal of the American Statistical Association,1975,70(350):386-393.

[11]Brorsen B W,Ryongyang S.Maximum likelihood estimates of symmetric stable distribution parameters[J].Communications on Statisties－Simulation,1990,19(4):1459-1464.

[12]石屹然,趙曉暉,單澤彪,等.α和高斯噪聲背景下線性極化陣列波達(dá)方向和極化參數(shù)聯(lián)合估計(jì)的FLOCC-SMN方法[J].吉林大學(xué)學(xué)報(bào):工學(xué)版,2016,46(4):1297-1303.Shi Yi-ran,Zhao Xiao-hui,Shan Ze-biao,et al.FLOCC-SMN method for DOA polarization parameters estimation of linear polarization array inαand Gaussian noise[J].Journal of Jilin University(Engineering and Technology Edition),2016,46(4):1297-1303.

[13]應(yīng)文威,歐勇恒,蔣宇中,等.新型自適應(yīng)非高斯接收機(jī)設(shè)計(jì)[J].吉林大學(xué)學(xué)報(bào):工學(xué)版,2013,43(6):1685-1689.Ying Wen-wei,Ou Yong-heng,Jiang Yu-zhong,et al.New adaptive receiver for channels with non-Gaussian noise[J].Journal of Jilin University(Engineering and Technology Edition),2013,43(6):1685-1689.

[14]劉文紅.脈沖噪聲下時(shí)間延遲估計(jì)方法及應(yīng)用的研究[D].大連:大連理工大學(xué)信息與通信工程學(xué)院,2007 Liu Wen-hong.Study on time delay estimation methods and applications in the presence of impulsive noise[D].Dalian:School of Information and Communication Engineering,Dalian University of Technology,2007.