對稱α穩定分布特性及其參數估計方法研究*

王 余 王本猛 衛紅凱 彭 成

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

對稱α穩定分布特性及其參數估計方法研究*

王 余1王本猛2衛紅凱3彭 成3

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

對稱α穩定分布能較好地描述混響等非高斯對稱概率密度序列特征，首先探討并給出序列非高斯性、對稱性判定方法，在此基礎上，研究并比較分數低階矩法和對數矩法兩種對稱α穩定分布參數估計方法，仿真及實測數據處理結果表明：對數矩法估計性能優于分數低階矩法。

對稱α穩定分布; 混響; 概率密度分布

Class Number TN911.7

1 引言

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

α穩定分布是廣義的高斯分布，具有四個參數，通過調整參數值可實現對穩定分布PDF拖尾厚度的控制，可以對PDF為單峰鐘形的隨機變量進行靈活的描述。若序列概率密度函數對稱，則可簡化α穩定分布的參數估計。本文首先探討了序列非高斯性、對稱性判定方法，并以仿真數據進行驗證；在此基礎上，研究了分數低階矩法和對數矩法兩種對稱α穩定分布參數估計方法，仿真及實測數據處理結果表明：對數矩法估計性能優于分數低階矩法。

2 α穩定分布

2.1 定義

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

若隨機變量X的特征函數具有如下形式：

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

(1)

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

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

1)α為特征指數，它決定了穩定分布脈沖特性的程度。α值越小，所對應穩定分布的拖尾越厚，表明穩定分布隨機變量X在遠離中心位置取值的概率越大，即脈沖特性越顯著。反之，則脈沖特性減弱。當α=2時，α穩定分布為高斯分布，即高斯分布是α穩定分布的一個特例。

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

3)β為偏斜指數，它決定了α穩定分布的偏斜程度。β=0時，α穩定分布是對稱的，稱為對稱α穩定分布，記為SαS(Symmtric alpha-stable distribution)。經典高斯分布(即正態分布)和柯西分布都屬于SαS分布。β>0和β<0分別對應于右偏斜和左偏斜。當α=1，β=0時，α穩定分布為柯西分布，即柯西分布是α穩定分布的一個特例。

2.2 α穩定分布的PDF

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

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

(2)

特別地，標準α穩定分布隨機變量的PDF為

(3)

可以證明，存在如下關系：

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

(4)

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

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

圖1～4給出了不同參數時的α穩定分布的PDF曲線。

3 序列非高斯性及對稱性確定方法

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

3.1 非高斯性判斷

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

(5)

(6)

3.2 對稱性判斷

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

對于方法1，取α=0.7，γ=1，μ=0，β分別為0、-1、1，得到10000個αS分布樣本序列點，繪制出相應的統計概率密度曲線圖，如圖9所示。

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

對于方法2，取α=0.7，γ=1，μ=0，β分別取0、±0.1、±0.8、±1，得到10000個αS分布樣本序列點，對其中的正負值個數進行統計，將統計結果記錄在表1中。

表1 樣本序列的正負值個數統計

從表1不難看出，β取不同值時，樣本序列的正負值統計個數是不相等的。特別的，當β=0時，即滿足對稱分布時，正負值個數基本相等。而當β取值越接近1時，正樣本值個數就越多，負樣本值個數就越少；當β取值越接近-1時，負樣本值個數就越多，正樣本值個數就越少。仿真結果與理論分析規律符合，因此，可依據統計出的正負樣本值個數是否近似相等判斷分布是否對稱，而不需事先知道β的具體取值。

4 α穩定分布的參數估計

α穩定分布的特性由其特征函數的4個參數α，β，γ和μ確定。而絕大多數水聲信號是對稱分布的，不失一般性，我們只討論參數α和γ的估計。

4.1 分數低階矩法

設X是一個實的SαS隨機變量，則其分數低階矩存在，如下式所示:

(5)

(6)

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

(7)

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

(8)

(9)

(10)

(11)

分數低階矩法需預先知道α的先驗信息，為了克服這一不足，引入更為一般的參數估計方法——對數矩法。

4.2 對數矩法

(12)

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

(13)

(14)

化簡后得到

(15)

其中，Ce=0.57721566…是Euler常數，且有

(16)

(17)

(18)

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

(19)

(20)

其中，N是樣本總數，Yi是獨立同分布的觀測值。

5 仿真/試驗數據估計

本節分別以仿真和試驗數據對上節中的兩種SαS參數估計方法進行比較。

5.1 仿真數據估計

5.1.1 分數低階矩法

1) 取α值為：0.2、0.8、1.6，p值為0.2。分別重復1000次產生10000個樣本點的SαS分布隨機變量序列，并對估計結果取平均值，則參數α和σ的估計結果如表2所示。

表2 α取不同值時分數低階矩法的估計結果

2) 取p值為：0.2、0.4、0.8，α值為0.6。分別重復1000次產生10000個樣本點的SαS分布隨機變量序列，并對估計結果取平均值，則參數α和σ的分數低階矩法估計結果如表3所示。

表3 p取不同值時分數低階矩法的估計結果

3) 樣本序列點數取不同值也會對估計結果準確性帶來影響。重復1000次產生α=0.6，σ=1的SαS分布序列，樣本序列點分別取1000、5000、10000，且p值取0.1。估計參數α和σ的值，并分別對不同情況下的估計結果求取平均值，分數低階矩法估計結果如表4所示。

觀察表2～4可知，分數低階矩法的估計結果比較準確，并且參數α的取值越小，其估計結果就越準確。p取值越小，其估計的結果也愈準確；然而當p的取值超過α時，估計的結果將會出現較大的誤差，因此為了得到準確的估計結果，盡量選擇較小的p值；樣本序列點數越大，對參數α和σ的估計越準確，估計的穩定性也越好。在通常情況下，樣本點數取為5000點就可滿足估計要求。

表4 樣本點數N取不同值時分數低階矩法的估計結果

5.1.2 對數矩法

1) 考察不同α取值下對數矩法性能。α值分別取：0.2、0.8、1.6，β=0，σ=1，μ=0。分別重復1000次產生10000個樣本點的SαS分布隨機變量序列，參數α和σ的對數矩法估計結果如表5所示。

表5 α取不同值時對數矩法的估計結果

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

表6 樣本點數N取不同值時對數矩法的估計結果

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

5.2 試驗數據估計

對數矩法不需先驗信息，且估計結果準確，以對數矩法對實測混響數據進行特征參量估計。由于特征指數α是影響混響的主要因素[12]，因此只對試驗數據的特征指數α進行估計。

取[0411三亞海試]和[0601三亞海試]不同基元級數據進行參數α估計，得到的估計結果如表7所示。原始數據均經相關中頻的512階FIR帶通濾波預處理。

表7 海洋實測數據的α估計值

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

6 結語

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

[1] 邱天爽,張旭秀,李小兵,等.統計信號處理：非高斯信號處理及其應用[M].北京：電子工業出版社,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日

王余,男,助理研究員,研究方向:海上水下戰略預警等。

TN911.7

10.3969/j.issn.1672-9730.2017.04.029