对称α稳定分布特性及其参数估计方法研究*

王 余 王本猛 卫红凯 彭 成

(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