α穩定分布特征指數估計算法

關濟實,石要武,邱建文,單澤彪,史紅偉

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

0 引 言

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

1 α穩定分布的定義和性質

1.1 α穩定分布的定義

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

定義 若一個隨機變量X的特征函數可以表示為:

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

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

1.2 α穩定分布的性質

現有文獻對α穩定分布的性質做了詳細的研究,本文不再一一說明,僅介紹本文算法相關的兩條性質,即疊加性質和p范數性質[14]。

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

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

2 基于組合分布和分數低階矩的系數估計方法

2.1 特征指數估計算法原理

根據α分布的性質1,兩個獨立的同特征指數的穩定分布變量相加時,其和變量同樣服從α分布,并且其分散系數仍為α,當多個同分布的α信號相加時,這一性質可以以如下形式表述:

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

對照性質1,此推理很容易得證。

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

綜合以上兩性質和推理,多個同分布變量相加時有:

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

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

由式(13)可知,當有多個同分布α穩定分布隨機變量時,可以將多個隨機變量組合,利用組合變量與原變量之間p范數關系,得到特征指數α的估計。

2.2 α穩定分布參數估計方法

2.2.1 參數估計方法

式(13)表明多個同分布的α信號相加時,其統計量E[|X|p]1/p之間的比只與參與相加的變量個數和特征指數α有關。利用這一特征,可以從實際采樣信號中估計得到特征指數α。

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

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

由以上定義可知,Y和Z是相互獨立的隨機變量。

算法的執行步驟如下:

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

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

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

(4)對X和X′分別計算統計量T和T′。

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

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

2.2.2 估計方法的無偏性與一致性分析

式(12)由α穩定分布的性質導出,因此它是嚴格成立的,但在實際運算中,無法直接得到變量的p階矩,只能以式(17)得到的估計,式(17)中統計量T′的數學期望為:

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

α與D之間是冪函數關系,二者之間是一一對應的,因此如果算法中對D的估計是無偏的,則對α的估計也將是無偏的,算法中D的估計值為:

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

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

3 數值仿真實驗

實驗1 不同α值的特征指數估計實驗

實驗目的:驗證算法的有效性。

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

估計結果與生成數據時所用α之間關系如圖1所示。

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

圖1 基本的α估計實驗結果Fig.1 Estimation result ofαwith basic experiment

每個α取值進行100次蒙特卡羅實驗,計算均方誤差,圖2為均方誤差圖。

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

可以看到,隨著α變大,均方誤差逐漸變大,當α＝2時,均方誤差僅為－12 dB,具有較高的估計精度。

實驗2 分割段數對算法的影響

實驗目的:驗證分割段數對于算法精度的影響。

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

圖3為分段數變化時特征指數的估計結果和均方誤差。

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

實驗3 分散系數對算法的影響實驗

實驗目的:驗證不同分散系數γ時的算法性能。

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

圖4為γ分段數變化時估計誤差結果。

圖3 分段數對α估計結果的影響Fig.3 Influence of section number on estimation ofα

圖4 分散系數對α估計結果的影響Fig.4 Influence ofγon estimation ofα

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

實驗4 高斯噪聲與SαS噪聲共存時的估計效果

實驗目的:驗證高斯信號與SαS信號共存時算法對于α的估計效果。

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

圖5為高斯噪聲摻雜時α的估計結果。

圖5 高斯噪聲摻雜時算法對α的估計結果Fig.5 Estimation ofαwhen Gaussian noise added

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

4 結束語

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

[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]孫永梅.α穩定分布參數估計與譜分析理論及應用研究[D].大連:大連理工大學信息與通信工程學院,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]李森.穩定分布噪聲下通信信號處理新算法及性能分析[D].大連:大連理工大學信息與通信工程學院,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]何勁.α穩定分布噪聲背景下陣列信號處理方法研究[D].南京:南京理工大學電子工程與光電技術學院,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]石屹然,趙曉暉,單澤彪,等.α和高斯噪聲背景下線性極化陣列波達方向和極化參數聯合估計的FLOCC-SMN方法[J].吉林大學學報:工學版,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]應文威,歐勇恒,蔣宇中,等.新型自適應非高斯接收機設計[J].吉林大學學報:工學版,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]劉文紅.脈沖噪聲下時間延遲估計方法及應用的研究[D].大連:大連理工大學信息與通信工程學院,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.