對(duì)稱α穩(wěn)定分布特性及其參數(shù)估計(jì)方法研究*

王 余 王本猛 衛(wèi)紅凱 彭 成

(1.31001部隊(duì) 北京 100080)(2.92001部隊(duì) 青島 266000)(3.海軍工程大學(xué) 武漢 430033)

對(duì)稱α穩(wěn)定分布特性及其參數(shù)估計(jì)方法研究*

王 余1王本猛2衛(wèi)紅凱3彭 成3

(1.31001部隊(duì) 北京 100080)(2.92001部隊(duì) 青島 266000)(3.海軍工程大學(xué) 武漢 430033)

對(duì)稱α穩(wěn)定分布能較好地描述混響等非高斯對(duì)稱概率密度序列特征，首先探討并給出序列非高斯性、對(duì)稱性判定方法，在此基礎(chǔ)上，研究并比較分?jǐn)?shù)低階矩法和對(duì)數(shù)矩法兩種對(duì)稱α穩(wěn)定分布參數(shù)估計(jì)方法，仿真及實(shí)測(cè)數(shù)據(jù)處理結(jié)果表明：對(duì)數(shù)矩法估計(jì)性能優(yōu)于分?jǐn)?shù)低階矩法。

對(duì)稱α穩(wěn)定分布; 混響; 概率密度分布

Class Number TN911.7

1 引言

混響作為一種非白非高斯背景干擾，基于傳統(tǒng)高斯分布假設(shè)抑制混響，經(jīng)常會(huì)由于模型與干擾背景噪聲不能很好匹配而導(dǎo)致所設(shè)計(jì)的信號(hào)處理器顯著退化。因此，放棄傳統(tǒng)高斯假設(shè)而采用非高斯假設(shè)在水聲信號(hào)處理業(yè)界已成為趨勢(shì)。近年來，能夠描述非高斯分布模型的α穩(wěn)定分布(Alpha Stable Distribution,αS)[1～9]，因其具有統(tǒng)計(jì)分布的穩(wěn)定性和概率密度函數(shù)(Probability Density Function,PDF)的代數(shù)拖尾特點(diǎn)，在信號(hào)處理領(lǐng)域受到了廣泛的重視，己成為常用的沖激噪聲數(shù)學(xué)模型。

α穩(wěn)定分布是廣義的高斯分布，具有四個(gè)參數(shù)，通過調(diào)整參數(shù)值可實(shí)現(xiàn)對(duì)穩(wěn)定分布PDF拖尾厚度的控制，可以對(duì)PDF為單峰鐘形的隨機(jī)變量進(jìn)行靈活的描述。若序列概率密度函數(shù)對(duì)稱，則可簡(jiǎn)化α穩(wěn)定分布的參數(shù)估計(jì)。本文首先探討了序列非高斯性、對(duì)稱性判定方法，并以仿真數(shù)據(jù)進(jìn)行驗(yàn)證；在此基礎(chǔ)上，研究了分?jǐn)?shù)低階矩法和對(duì)數(shù)矩法兩種對(duì)稱α穩(wěn)定分布參數(shù)估計(jì)方法，仿真及實(shí)測(cè)數(shù)據(jù)處理結(jié)果表明：對(duì)數(shù)矩法估計(jì)性能優(yōu)于分?jǐn)?shù)低階矩法。

2 α穩(wěn)定分布

2.1 定義

一般來說，可以從穩(wěn)定性(Stability Property)、吸收域(Domain of Attraction)和特征函數(shù)(Characteristic Function)三個(gè)方面對(duì)α穩(wěn)定分布進(jìn)行定義[10]。其中，基于特征函數(shù)的α穩(wěn)定分布應(yīng)用最為廣泛。

若隨機(jī)變量X的特征函數(shù)具有如下形式：

φ(t)=E[exp(itX)]

(1)

則稱隨機(jī)變量X服從α穩(wěn)定分布。其中，β∈[-1,1]，μ∈(-∞,∞)。

α穩(wěn)定分布由參數(shù)α，β，γ，μ唯一確定，因此，將服從α穩(wěn)定分布的隨機(jī)變量X記為X∽S(α,β,γ,μ)。參數(shù)α，β，γ，μ具有如下物理意義：

1)α為特征指數(shù)，它決定了穩(wěn)定分布脈沖特性的程度。α值越小，所對(duì)應(yīng)穩(wěn)定分布的拖尾越厚，表明穩(wěn)定分布隨機(jī)變量X在遠(yuǎn)離中心位置取值的概率越大，即脈沖特性越顯著。反之，則脈沖特性減弱。當(dāng)α=2時(shí)，α穩(wěn)定分布為高斯分布，即高斯分布是α穩(wěn)定分布的一個(gè)特例。

2)γ為離差，也稱為尺度參數(shù)(Scale Parameter)，它決定穩(wěn)定分布隨機(jī)變量偏離其均值或中值的程度。γ=ζα，稱ζ為標(biāo)準(zhǔn)離差。γ類似于高斯分布中的方差，同樣，ζ類似于高斯分布中的標(biāo)準(zhǔn)偏差σ(但不相等)。在高斯情況下，γ的數(shù)值是方差σ2的一半，即有σ2=2γ。

3)β為偏斜指數(shù)，它決定了α穩(wěn)定分布的偏斜程度。β=0時(shí)，α穩(wěn)定分布是對(duì)稱的，稱為對(duì)稱α穩(wěn)定分布，記為SαS(Symmtric alpha-stable distribution)。經(jīng)典高斯分布(即正態(tài)分布)和柯西分布都屬于SαS分布。β>0和β<0分別對(duì)應(yīng)于右偏斜和左偏斜。當(dāng)α=1，β=0時(shí)，α穩(wěn)定分布為柯西分布，即柯西分布是α穩(wěn)定分布的一個(gè)特例。

2.2 α穩(wěn)定分布的PDF

α穩(wěn)定分布特征函數(shù)是其PDF的傅立葉變換，因此可通過對(duì)α穩(wěn)定分布隨機(jī)變量的特征函數(shù)求取傅立葉逆變換得到其PDF，如式(2)所示：

cos[(x-μ)t+γtαβω(t,α)]dt

(2)

特別地，標(biāo)準(zhǔn)α穩(wěn)定分布隨機(jī)變量的PDF為

(3)

可以證明，存在如下關(guān)系：

f(x;α,β,γ,μ)=f(-x;α,-β,γ,-μ)

(4)

即α穩(wěn)定分布的PDF相對(duì)于隨機(jī)變量x、偏斜指數(shù)β和位置參數(shù)μ是偶對(duì)稱的。另一方面，還可以證明f(x;α,γ,β,μ)是有界的，且任意階導(dǎo)數(shù)存在。

除高斯分布(Gaussian：α=2，β=0，σ2=2γ)、柯西分布(Cauchy：α=1，β=0)、列維分布(Levy：α=1/2，β=1)以外，一般穩(wěn)定分布的PDF均沒有封閉的表達(dá)式，需要通過數(shù)值積分方式計(jì)算其PDF。

圖1～4給出了不同參數(shù)時(shí)的α穩(wěn)定分布的PDF曲線。

3 序列非高斯性及對(duì)稱性確定方法

在一般應(yīng)用情況中，首先要確定所得序列的分布類型，才能基于相應(yīng)的分布性質(zhì)展開研究。在判定序列分布是否滿足非高斯特性基礎(chǔ)上，再進(jìn)一步假設(shè)該序列滿足αS分布。比如，對(duì)于自適應(yīng)濾波[11～12]等應(yīng)用，需在判斷信號(hào)是否滿足非高斯分布基礎(chǔ)上進(jìn)行。

3.1 非高斯性判斷

如果xk,k=1,2,…,N是一串滿足αS分布的樣本序列值，對(duì)于任何一個(gè)1≤n≤N范圍內(nèi)的n，其前n個(gè)樣本值的方差和均值分別為

(5)

(6)

3.2 對(duì)稱性判斷

αS分布可分為斜分布和對(duì)稱分布兩類，不同類分布下的參數(shù)估計(jì)方法不盡相同。在αS分布的特征函數(shù)中，偏斜指數(shù)β決定分布的偏斜程度。SαS分布是β=0情況時(shí)的αS分布，此時(shí)對(duì)除了α外的其他參數(shù)的估計(jì)將得到簡(jiǎn)化。所以在對(duì)αS分布的各參數(shù)進(jìn)行估計(jì)之前，預(yù)先判斷αS分布的對(duì)稱性是十分重要的。這里給出兩種對(duì)稱性判斷方法：方法1是直接繪制出樣本序列的統(tǒng)計(jì)PDF曲線圖，對(duì)其對(duì)稱性進(jìn)行直觀判斷；方法2是如果αS分布滿足對(duì)稱性，則其正負(fù)樣本序列值個(gè)數(shù)近似相等，此時(shí)可根據(jù)正負(fù)樣本序列值個(gè)數(shù)是否近似相等來判斷樣本的對(duì)稱性。

對(duì)于方法1，取α=0.7，γ=1，μ=0，β分別為0、-1、1，得到10000個(gè)αS分布樣本序列點(diǎn)，繪制出相應(yīng)的統(tǒng)計(jì)概率密度曲線圖，如圖9所示。

從圖9可看出，β取值為0時(shí)，PDF曲線具有對(duì)稱性；β取值為1時(shí)，PDF曲線往右偏；β取值為-1時(shí)，PDF曲線往左偏。所以可通過觀察樣本序列的統(tǒng)計(jì)PDF曲線確定是否滿足對(duì)稱分布。

對(duì)于方法2，取α=0.7，γ=1，μ=0，β分別取0、±0.1、±0.8、±1，得到10000個(gè)αS分布樣本序列點(diǎn)，對(duì)其中的正負(fù)值個(gè)數(shù)進(jìn)行統(tǒng)計(jì)，將統(tǒng)計(jì)結(jié)果記錄在表1中。

表1 樣本序列的正負(fù)值個(gè)數(shù)統(tǒng)計(jì)

從表1不難看出，β取不同值時(shí)，樣本序列的正負(fù)值統(tǒng)計(jì)個(gè)數(shù)是不相等的。特別的，當(dāng)β=0時(shí)，即滿足對(duì)稱分布時(shí)，正負(fù)值個(gè)數(shù)基本相等。而當(dāng)β取值越接近1時(shí)，正樣本值個(gè)數(shù)就越多，負(fù)樣本值個(gè)數(shù)就越少；當(dāng)β取值越接近-1時(shí)，負(fù)樣本值個(gè)數(shù)就越多，正樣本值個(gè)數(shù)就越少。仿真結(jié)果與理論分析規(guī)律符合，因此，可依據(jù)統(tǒng)計(jì)出的正負(fù)樣本值個(gè)數(shù)是否近似相等判斷分布是否對(duì)稱，而不需事先知道β的具體取值。

4 α穩(wěn)定分布的參數(shù)估計(jì)

α穩(wěn)定分布的特性由其特征函數(shù)的4個(gè)參數(shù)α，β，γ和μ確定。而絕大多數(shù)水聲信號(hào)是對(duì)稱分布的，不失一般性，我們只討論參數(shù)α和γ的估計(jì)。

4.1 分?jǐn)?shù)低階矩法

設(shè)X是一個(gè)實(shí)的SαS隨機(jī)變量，則其分?jǐn)?shù)低階矩存在，如下式所示:

(5)

(6)

其中，α是特征指數(shù)(0<α<2)，γ是分散系數(shù)，γ=σα，Γ(·)為伽馬函數(shù)，定義如下：

(7)

由式(5)～式(7)，可得到參數(shù)α的估計(jì)表達(dá)式如下所示：

(8)

(9)

(10)

(11)

分?jǐn)?shù)低階矩法需預(yù)先知道α的先驗(yàn)信息，為了克服這一不足，引入更為一般的參數(shù)估計(jì)方法——對(duì)數(shù)矩法。

4.2 對(duì)數(shù)矩法

(12)

其中，E(epY)是Y的矩生成函數(shù)，如下所示：

(13)

(14)

化簡(jiǎn)后得到

(15)

其中，Ce=0.57721566…是Euler常數(shù)，且有

(16)

(17)

(18)

則Y的期望和方差可通過以下兩式估計(jì)：

(19)

(20)

其中，N是樣本總數(shù)，Yi是獨(dú)立同分布的觀測(cè)值。

5 仿真/試驗(yàn)數(shù)據(jù)估計(jì)

本節(jié)分別以仿真和試驗(yàn)數(shù)據(jù)對(duì)上節(jié)中的兩種SαS參數(shù)估計(jì)方法進(jìn)行比較。

5.1 仿真數(shù)據(jù)估計(jì)

5.1.1 分?jǐn)?shù)低階矩法

1) 取α值為：0.2、0.8、1.6，p值為0.2。分別重復(fù)1000次產(chǎn)生10000個(gè)樣本點(diǎn)的SαS分布隨機(jī)變量序列，并對(duì)估計(jì)結(jié)果取平均值，則參數(shù)α和σ的估計(jì)結(jié)果如表2所示。

表2 α取不同值時(shí)分?jǐn)?shù)低階矩法的估計(jì)結(jié)果

2) 取p值為：0.2、0.4、0.8，α值為0.6。分別重復(fù)1000次產(chǎn)生10000個(gè)樣本點(diǎn)的SαS分布隨機(jī)變量序列，并對(duì)估計(jì)結(jié)果取平均值，則參數(shù)α和σ的分?jǐn)?shù)低階矩法估計(jì)結(jié)果如表3所示。

表3 p取不同值時(shí)分?jǐn)?shù)低階矩法的估計(jì)結(jié)果

3) 樣本序列點(diǎn)數(shù)取不同值也會(huì)對(duì)估計(jì)結(jié)果準(zhǔn)確性帶來影響。重復(fù)1000次產(chǎn)生α=0.6，σ=1的SαS分布序列，樣本序列點(diǎn)分別取1000、5000、10000，且p值取0.1。估計(jì)參數(shù)α和σ的值，并分別對(duì)不同情況下的估計(jì)結(jié)果求取平均值，分?jǐn)?shù)低階矩法估計(jì)結(jié)果如表4所示。

觀察表2～4可知，分?jǐn)?shù)低階矩法的估計(jì)結(jié)果比較準(zhǔn)確，并且參數(shù)α的取值越小，其估計(jì)結(jié)果就越準(zhǔn)確。p取值越小，其估計(jì)的結(jié)果也愈準(zhǔn)確；然而當(dāng)p的取值超過α?xí)r，估計(jì)的結(jié)果將會(huì)出現(xiàn)較大的誤差，因此為了得到準(zhǔn)確的估計(jì)結(jié)果，盡量選擇較小的p值；樣本序列點(diǎn)數(shù)越大，對(duì)參數(shù)α和σ的估計(jì)越準(zhǔn)確，估計(jì)的穩(wěn)定性也越好。在通常情況下，樣本點(diǎn)數(shù)取為5000點(diǎn)就可滿足估計(jì)要求。

表4 樣本點(diǎn)數(shù)N取不同值時(shí)分?jǐn)?shù)低階矩法的估計(jì)結(jié)果

5.1.2 對(duì)數(shù)矩法

1) 考察不同α取值下對(duì)數(shù)矩法性能。α值分別?。?.2、0.8、1.6，β=0，σ=1，μ=0。分別重復(fù)1000次產(chǎn)生10000個(gè)樣本點(diǎn)的SαS分布隨機(jī)變量序列，參數(shù)α和σ的對(duì)數(shù)矩法估計(jì)結(jié)果如表5所示。

表5 α取不同值時(shí)對(duì)數(shù)矩法的估計(jì)結(jié)果

2) 考察不同樣本點(diǎn)數(shù)N情況下對(duì)數(shù)矩法估計(jì)性能。取α=1.2，β=0，σ=1，μ=0，樣本點(diǎn)數(shù)N分別取1000，5000，10000，分別1000次產(chǎn)生不同點(diǎn)數(shù)的SαS分布序列，參數(shù)α和σ的對(duì)數(shù)矩法估計(jì)結(jié)果如表6所示。

表6 樣本點(diǎn)數(shù)N取不同值時(shí)對(duì)數(shù)矩法的估計(jì)結(jié)果

比較表2和表5、表4和表6可知，分?jǐn)?shù)低階矩法與對(duì)數(shù)矩法的參數(shù)估計(jì)性能相似，都能較有效地估計(jì)出參數(shù)值，并且參數(shù)α的取值越小，估計(jì)出的結(jié)果就越準(zhǔn)確，樣本點(diǎn)數(shù)越多，估計(jì)出的結(jié)果也更加準(zhǔn)確，兩種參數(shù)估計(jì)方法中樣本點(diǎn)數(shù)N取5000時(shí)估計(jì)精度已足夠滿足實(shí)際需求。由于分?jǐn)?shù)低階矩法需要通過解sinc方程才能得到估計(jì)結(jié)果，而對(duì)數(shù)矩法既不需要解sinc方程也不需要額外確定參數(shù)p的取值(即不需要α的先驗(yàn)信息)，同時(shí)參數(shù)估計(jì)公式也更簡(jiǎn)單，且其估計(jì)結(jié)果更準(zhǔn)確，因此對(duì)數(shù)矩法的魯棒性更強(qiáng)。

5.2 試驗(yàn)數(shù)據(jù)估計(jì)

對(duì)數(shù)矩法不需先驗(yàn)信息，且估計(jì)結(jié)果準(zhǔn)確，以對(duì)數(shù)矩法對(duì)實(shí)測(cè)混響數(shù)據(jù)進(jìn)行特征參量估計(jì)。由于特征指數(shù)α是影響混響的主要因素[12]，因此只對(duì)試驗(yàn)數(shù)據(jù)的特征指數(shù)α進(jìn)行估計(jì)。

取[0411三亞海試]和[0601三亞海試]不同基元級(jí)數(shù)據(jù)進(jìn)行參數(shù)α估計(jì)，得到的估計(jì)結(jié)果如表7所示。原始數(shù)據(jù)均經(jīng)相關(guān)中頻的512階FIR帶通濾波預(yù)處理。

表7 海洋實(shí)測(cè)數(shù)據(jù)的α估計(jì)值

水下混響背景噪聲滿足SαS分布模型，隨著水下環(huán)境噪聲強(qiáng)度的增加，背景噪聲的非高斯脈沖性降低，相應(yīng)的SαS分布模型的特征指數(shù)α也將增大。由文獻(xiàn)[8]可知，大部分海洋混響噪聲的特征指數(shù)α取值集中在1.6～2之間。結(jié)合表7對(duì)三亞海洋測(cè)試數(shù)據(jù)特征指數(shù)α的估計(jì)結(jié)果可以知道，實(shí)際海洋混響噪聲的特征指數(shù)α取值基本在1.65～1.9之間。

6 結(jié)語

本文首先對(duì)序列的非高斯性與對(duì)稱性的判斷進(jìn)行了探討，并結(jié)合仿真試驗(yàn)給出了判定方法；在此基礎(chǔ)上，研究并比較了兩種SαS模型參數(shù)估計(jì)方法，即分?jǐn)?shù)低階矩法和對(duì)數(shù)矩法，仿真和實(shí)測(cè)數(shù)據(jù)估計(jì)結(jié)果表明：兩種方法均能較準(zhǔn)確地估計(jì)SαS模型參數(shù)值，但分?jǐn)?shù)低階矩法需事先知道特征指數(shù)α的先驗(yàn)值，而對(duì)數(shù)矩法不需先驗(yàn)信息，且估計(jì)方法簡(jiǎn)單，普適性更好。

[1] 邱天爽,張旭秀,李小兵,等.統(tǒng)計(jì)信號(hào)處理：非高斯信號(hào)處理及其應(yīng)用[M].北京：電子工業(yè)出版社,2004,20-100.[2] Shao M, Nikias C L. Signal processing with fractional lower order moments: stable processes and their applications [J]. Proceedings of the IEEE, 1993, 81(7): 986-1010.

[3] Rupi M, Tsakalides P, Nikias C L, et al. Robust constant modulus array based on fractional lower order statistics [C]//IEEE International Conference on Acoustics, Speech, and Signal Processing,1999,5:2945-2948.

[4] Gonzalez J G, Griffith D W, Arce G R. Zero-order statistics: a signal processing framework for very impulsive processes[J]. Processings of the IEEE Signal Processing Workshop on Higher-Order Statistics,1997:254-258.

[5] Kuruoglu E E. Analytical representation for positive α-stable densities [C]//ICASSP 2003,2003,6:729-732.

[6] Tsihrintzis G A, Nikias C L. Data-adaptive algorithms for signal detection in sub-Gussian impulsive interference [J]. IEEE Transactions on Signal Processing,1997,45(7): 1873-1878.

[7] Bondenschatz J S, Nikias C L. Symmetric alpha-stable filter theory [C]//IEEE Transactions on Signal Processing,1997,45(9):2301-2306.

[8] Xinyu Ma, Nikias C L. Joint estimation of time delay and frequency delay in impulsive noise using fractional lower order statistics[J]. IEEE Transactions on Signal Processing, 1996,44(11):2669-2687.

[9] Shao M, Nikias C L. Detection and adaptive estimation of stable processes with fractional lower-order moments[C]//IEEE Sixth SP Workshop on Statistical Signal and Array Processing,1992,94-97.

[10] Manolakis D G, Ingle V K, Kogon S M. Statistical and adaptive signal processing[M]. Boston: McGraw-Hill,2000,40-80.

[11] M.Belge, E.L.Miller. A sliding window RLS-like adaptive algorithm for filtering alpha-stable noise[J]. IEEE Signal Processing Letters,2000,7(4):86-89.

[12] Shao M, Nikias C L. Signal detection in impulsive noise based on stable distributions [C]//1993 Conference Record of The Twenty-Seventh Asilomar Conference on Signals, Systems and Computers, 1993, 1: 218-222.

Symmetrical Alpha Stable Distribution and Its Parameter Estimation Method

WANG Yu1WANG Benmeng2WEI Hongkai3PENG Cheng3

(1. No. 31001 Troops of PLA, Beijing 100080)(2. No. 92001 Troops of PLA, Qingdao 266000)(3. Naval University of Engineering, Wuhan 430033)

Symmetrical alpha stable distribution can describe sequence with symmetry and non-Gaussian probability density such as reverberation. The paper first discusses and gives the method of judging symmetry and non-Gaussian property of probability density. Then based on this, two parameter estimation methods for symmetrical alpha stable distribution that are fractional lower-order moment method and logarithmic moment method are studied and compared. Simulation and experimental results show that the logarithmic moment method is preferable to fractional lower-order moment method.

symmetry alpha stable distribution, reverberation, probability density distribution

2016年10月11日,

2016年11月29日

王余,男,助理研究員,研究方向:海上水下戰(zhàn)略預(yù)警等。

TN911.7

10.3969/j.issn.1672-9730.2017.04.029