在传统的声呐信号处理中，高斯模型占据着主导地位。随着现代主动声呐体制向低频、大功率、大孔径方向发展，混响对声呐性能的影响日趋严重，如何有效地抑制混响干扰成为主动声呐信号处理的首要任务。大量实践表明，沿袭传统的高斯分布假设，往往不能取得理想的抗混响效果。因此，放弃高斯假设而采用非高斯分布假设在水声信号处理业界已成为一种趋势^[1-2]。近年来，α稳定分布作为一种非高斯分布模型，具有统计分布的稳定性和概率密度函数(probability density function，PDF)的代数拖尾特点，在信号处理领域受到了广泛的重视，己成为常用的冲激信号统计模型^[3-7]。作为唯一的一类构成独立同分布随机变量之和的极限分布^[4]，α稳定分布是广义的高斯分布，通过调整特征指数α的值实现对PDF拖尾厚度的控制，可以对PDF为单峰钟形随机信号统计特性进行灵活地描述，比高斯分布具有更广泛的适用性。而以混响为主的主动声呐背景干扰，其瞬时幅度样本是均值为零、正负样本量大致相等、上下包络线基本对称的随机信号^[1-2]，即其PDF恰好是以0为中心左右对称的单峰钟形曲线。因此，本文将原点位置对称α稳定分布(SαS, symmetric alpha-stable distribution)引入到混响非高斯PDF建模中，并使用仿真和实测数据对这一建模方法的可行性和正确性进行检验，最后通过大量海试混响数据建模结果分析，对混响数据的SαS分布规律进行归纳总结。

1 混响的特点与SαS分布简介 1.1 混响数据的特点

海洋混响是海洋信道中一切散射作用的总和。从瞬时幅度上看，它是紧跟在发射之后、总体由强至弱的一连串衰减(可能还有起伏震荡)波形。

图 1给出了一段典型的主动声呐阵元级接收数据波形，从中可以分析归纳混响数据的一般特点。

引发这段数据的发射信号是4 s脉宽的HFM脉冲。前4 s为发射脉冲直达声，称为直达声区。事实上，发射伊始即引发混响，亦即，直达声区中也有不完整发射脉冲引发的混响，但因其混叠于直达声中难以分离，且此处视为主动声呐近程盲区，一般不对其内的混响进行研究。直达声结束后，数据幅度从与直达声等高度逐渐衰减，直至不再明显降低，大致如图中的4~12 s段所示，这一段是典型混响区域，称为混响限制区，该区域内，混响强度高于背景噪声，是声呐的主要背景干扰。大致从12 s向后，混响淹没于背景噪声之中，以海洋环境噪声为主的背景噪声成为主要的背景干扰，混响的影响可以忽略，而背景噪声是宽平稳的，故接收数据幅度不再明显变化，此段称为噪声限制区。根据声呐量程，噪声限制区可以延续很长，这里仅给出一部分。不难理解，研究混响数据的统计特性，要基于混响限制区内的数据进行。

从图 1中还直观可见，混响数据是0均值的、正负样本量大致相等、上下包络线基本对称，这意味着，可以使用0均值对称单钟型曲线来对混响数据的概率密度函数进行拟合建模。零均值高斯分布、混合高斯分布，和下面要介绍的SαS分布，其PDF恰好就是这样一类曲线。

1.2 SαS分布简介

若随机变量X具有如下形式的特征函数^[3]：

$ \begin{array}{*{20}{c}} {\varphi (t) = E[\exp ({\rm{i}}tX)] = }\\ {\exp \left\{ {{\rm{i}}\mu t - \gamma |t{|^\alpha }[1 + {\rm{i}}\beta {\mathop{\rm sgn}} (t)\omega (t,\alpha )]} \right\}} \end{array} $

(1)

则称X服从α稳定分布，记为$X \sim S\left( {\alpha , \beta , \gamma , {\rm{ }}\mu } \right)$。其中，

$ \omega (t,\alpha ) = \left\{ {\begin{array}{*{20}{l}} {\tan \left( {\frac{{{\rm{ \mathsf{ π} }}\alpha }}{2}} \right),\alpha \ne 1}\\ { - \frac{2}{{\rm{ \mathsf{ π} }}}\ln |t|,\alpha = 1} \end{array},{\mathop{\rm sgn}} (t) = \left\{ {\begin{array}{*{20}{c}} {1,}&{t > 0}\\ {0,}&{t = 0}\\ { - 1,}&{t < 0} \end{array}} \right.} \right. $

式中：α为特征指数，表征稳定分布脉冲特性强弱，α∈(0, 2]，值越小，则所对应分布的拖尾越厚，非高斯脉冲特性越显著；γ为离差，亦称为尺度参数、分散系数，表征稳定分布随机变量偏离其均值或中值的程度，γ∈(0, ∞)，定义σ=γ^1/α为标准离差；β为偏斜指数，表这种分布称为原点征稳定分布的偏斜程度，β∈[－1, 1]，β=0、β>0和β<0时，分别对应稳定分布无偏斜(对称)、右偏斜和左偏斜情况；μ为位置参数，表征稳定分布的均值(α>1时)或中值(α≤1时)，μ∈(－∞, ∞)。

α稳定分布特征函数与其PDF是一对傅立叶变换^[3]，可通过对特征函数进行傅立叶变换得到PDF：

$ \begin{array}{*{20}{c}} {f(x;\alpha ,\beta ,\gamma ,\mu ) = }\\ {\frac{1}{{\rm{ \mathsf{ π} }}}\int_0^\infty {\exp } \left\{ { - \gamma {t^\alpha }} \right\}\cos \left[ {(x - \mu )t + \gamma {t^\alpha }\beta \omega (t,\alpha )} \right]{\rm{d}}t} \end{array} $

(2)

除几组特殊参数取值外，式(2)一般没有封闭的解析表达式，需要以数值积分方式计算^[8-9]。图 2给出了以这种方式计算绘出的不同参数取值时的α稳定分布的PDF曲线比较。从各图中参数变化时的PDF曲线对比，可以清晰地看到对应参数的物理意义体现。

一些常见分布可以看作是α稳定分布S(α, β, γ, μ)的特例^[3]，比如：当α=2、β=0时，S(2, 0, γ, μ)就是高斯分布(μ, 2γ)；当α=1、β=0时，S(1, 0, γ, μ)就是柯西分布(μ, γ)。

当偏斜指数β=0时，α稳定分布称为对称α稳定分布，记为SαS分布。不难发现，高斯分布、柯西分布都是SαS分布的特例。正常状态下的主动声呐背景干扰数据，包括混响数据，其样本都是无偏斜且0均值的，故可用位置参数μ=0的SαS分布来拟合其PDF，可把这种分布称为原点位置SαS分布，简记作Sα(γ)。

可以看到，S_α(γ)分布仅有2个参数：特征指数α和离差γ，就像常用的高斯分布一样简洁。值得指出，在对混响数据进行非高斯PDF建模时，堪用的模型，即使最简洁的二阶0均值混合高斯模型ZMGM(2)，也有3个参数^[1]。相对少的参数意味着相对快捷而精确的参数估计算法，在这一点上，S_α(γ)具有无可比拟的优势。

2 混响数据SαS分布建模方法

获取合格而足量的混响样本数据后，即可以原点位置SαS分布对其进行建模。而建模的核心任务是，如何构建S_α(γ)的2个未知参数α和γ估计器。

一般的α稳定分布S(α, β, γ, μ)参数估计问题，已经有很多文献研究^{[3, 8, 10]}。如前分析，混响数据实际上可以简化的α稳定分布，即S_α(γ) 分布进行PDF建模，故针对这一特点，本文从诸多方法中予以选取并简化改进，建立了如下2种方法：正负矩法和对数矩法。

2.1 正负矩法

设随机变量X~S_α(γ)，则有^[3]：

$ \begin{array}{l} {\rm{E}}\left( {|X{|^p}} \right){\rm{E}}\left( {|X{|^{ - p}}} \right) = C(p,\alpha )C( - p,\alpha ) = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{2\tan (p{\rm{ \mathsf{ π} }}/2)}}{{\alpha \sin (p{\rm{ \mathsf{ π} }}/\alpha )}},0 < p < \min (\alpha ,1) \end{array} $

(3)

式中：

$ {{M_{|X|}}(p) \buildrel \Delta \over = {\rm{E}}\left( {|X{|^p}} \right) = C(p,\alpha ){\gamma ^{p/\alpha }}} $

(4)

$ {C(p,\alpha ) = \frac{{{2^{p + 1}}\varGamma [(p + 1)/2]\varGamma ( - p/\alpha )}}{{\alpha \sqrt {\rm{ \mathsf{ π} }} \varGamma ( - p/2)}}} $

(5)

这里，Γ(·)为伽马函数，定义为^[8]：

$ \varGamma (x) \buildrel \Delta \over = \int_0^\infty {{t^{x - 1}}} {{\rm{e}}^{ - t}}{\rm{d}}t $

(6)

由式(3)~(5)，有：

$ \begin{array}{*{20}{c}} {{\mathop{\rm sinc}\nolimits} \left( {\frac{{p{\rm{ \mathsf{ π} }}}}{\alpha }} \right) = \frac{{\sin (p{\rm{ \mathsf{ π} }}/\alpha )}}{{(p{\rm{ \mathsf{ π} }}/\alpha )}} = }\\ {\frac{{2\tan (p{\rm{ \mathsf{ π} }}/2)}}{{p{\rm{ \mathsf{ π} }}{M_{|X|}}(p){M_{|X|}}( - p)}},0 < p < \min (\alpha ,1)} \end{array} $

(7)

通过样本统计得到$ \left| X \right|$ 的p阶和－p阶矩估计$ {\hat M_{\left| X \right|}}(p)$ 和${\hat M_{\left| X \right|}}( - p)$ 后，代入(7)即可解出特征指数估计值 $ \hat \alpha$ ，再代入式(4)即可解出离差估计值$ \hat \gamma$，完成S_α(γ)分布参数估计。不难发现，这种方法的起点是$ \left| X \right|$ 的正、负p阶矩估计，故称为正负矩法，简记为NM。

阶数p的设定是正负矩法的关键之一，但多数文献中仅给出一范围指示^{[3, 10]}，并无深入探讨。针对混响建模实际，本文采取如下做法：首先要根据经验确定待建模数据的特征指数下限值α_min(这一经验值的获取可通过下文所述的对数矩法建模总结)，确保p取值小于α_min或1；为提高估计精度，可设置多组阶数值p_i通过式(7)得到多组估值${\hat \alpha _i}$，进行平均得到最终估值。

2.2 对数矩法

设随机变量X~S_α(γ)，则由式(4)有：

$ {\rm{E}}\left( {|X{|^p}} \right) = {\rm{E}}\left( {{{\rm{e}}^{p\ln |X|}}} \right) = C(p,\alpha ){\gamma ^{p/\alpha }} $

(8)

引入负阶矩和分数阶矩概念后，式(8)在p=0处连续^[3]。定义随机变量Y=ln|X|，则Y的矩母函数为：

$ {\rm{E}}\left( {{{\rm{e}}^{pY}}} \right) \buildrel \Delta \over = \sum\limits_{k = 0}^\infty {\rm{E}} \left( {{Y^k}} \right)\frac{{{p^k}}}{{k!}} = C(p,\alpha ){\gamma ^{p/\alpha }} $

(9)

故Y的任意k阶矩都是有限的，且满足：

$ {\rm{E}}\left( {{Y^k}} \right) = {\left. {\frac{{{{\rm{d}}^k}}}{{{\rm{d}}{p^k}}}\left[ {C(p,\alpha ){\gamma ^{p/\alpha }}} \right]} \right|_{p = 0}} $

(10)

k=1时可得Y 的一阶原点矩(均值)为:

$ {\mu _Y} = {\rm{E}}(Y) = \left( {\frac{1}{\alpha } - 1} \right){c_e} + \frac{1}{\alpha }\ln \gamma $

(11)

式中c_e为Euler常数。同理可得Y的二阶中心距(方差)为：

$ \sigma _Y^2 = {\rm{E}}\left\{ {{{[Y - {\rm{E}}(Y)]}^2}} \right\} = \frac{{{{\rm{ \mathsf{ π} }}^2}}}{6}\left( {\frac{1}{{{\alpha ^2}}} + \frac{1}{2}} \right) $

(12)

混响建模实际应用中，可首先根据样本统计得到Y=ln|X|的均值和方差估计$\hat \mu$ _Y和$ \hat \sigma$ _Y²，然后把$ \hat \sigma$_Y²值代入式(12)解得特征指数估计值$ {\hat \alpha _i}$，最后再根据式(11)解得离差估计值 $\hat \gamma$ ，完成S_α(γ)分布参数估计。可以看到，这种方法的起点是 $\left| X \right| $ 的对数变换及其一、二阶矩估计，故称为对数矩法，简记为LM。

应当指出，式(11)与(12)中已无矩母函数自变量(亦即正负矩法的中阶数)p，这就是说，对数矩法不需要p初值设定，恰与正负矩法形成补充，二者估计结果可以相互印证。

2.3 仿真实例

首先使用PDF可设定的仿真数据检验S_α(γ) 分布对高斯和非高斯PDF的拟合能力。

图 3是一段标准正态分布随机序列的建模拟合情况。其中图 3(a)为波形图，图 3(b)为PDF曲线对比，虚线为使用设定的高斯分布参数[μ, σ²] =[0, 1]绘出的曲线，代表PDF真值；实线为使用对数矩法S_α(γ)分布参数估计[$ \hat \alpha , \hat \gamma$ ] =[2.0, 0.5]绘出的曲线，代表PDF拟合值。2条曲线几乎完全重合，这说明对于高斯数据，S_α(γ) 具有良好的PDF拟合性能。

图 4是一段非高斯随机序列的建模拟合情况。非高斯序列按参数[ε_I, σ_B², σ_I²] =[0.1, 1, 100]的混合高斯分布以复合抽样法^[1]生成。(a)为波形图，从中可观察到大量“毛刺”(尖峰脉冲)突出，故这是一种比较严重的非高斯情形。图 4(b)为PDF曲线对比，虚线为使用设定的混合高斯参数绘出的曲线，代表PDF真值；细实线为使用对数矩法S_α(γ)分布参数估计($ \hat \alpha , \hat \gamma$) =[1.4, 0.7]绘出的曲线，代表PDF拟合值。可以看到，这两条曲线几乎完全重合，这说明对于非高斯数据，S_α(γ) 同样具有良好的PDF拟合性能。为了比较，图 4(b)图中还绘出了以高斯分布拟合的PDF曲线，参数估计为[$ \hat \mu , \hat \sigma$ ²]=[0, 10.1]，如粗实线所示，可见，这条曲线已远远偏离PDF真值，说明高斯建模在这种情况下已不可用。

3 混响数据的SαS分布特性分析

本节以数组低频主动拖曳线列阵声呐(0411三亚海试)数据分析为例，验证以混响为主的主动声呐背景干扰数据(即图 1所示中的混响限制区内数据，简称为混响)SαS分布PDF建模方法的正确性，并初步归纳混响PDF的SαS分布规律。

3.1 建模方法验证

(0411三亚海试)发射3种主动脉冲信号：Ⅰ表示脉宽为432 ms、频率为750 Hz的CW脉冲，称为窄带短脉冲；Ⅱ表示同样脉宽、频带为650~850 Hz的LFM脉冲，称为宽带短脉冲；Ⅲ表示脉宽为1 296 ms、同样频带的LFM脉冲，称为宽带长脉冲。使用80元拖曳线列阵接收回波数据。进行PDF建模前数据已经过带通滤波→波束形成→移频→降采样→幅值归一化等处理^{[1, 11]}。

图 5~7分别给出了3段3种脉冲混响数据的波形图和SαS分布PDF建模拟合图。为便于比较，图 5~7(b)图中同时绘出了：①使用柱状图统计法得到的PDF曲线，图中标记为HS。在足够样本量下，HS曲线可视为真实分布。②使用正负矩法S_α(γ) 参数估计得到PDF拟合曲线，图中标记为NM。③使用对数矩法S_α(γ) 参数估计得到PDF拟合曲线，图中标记为LM。NM和LM方法得到的α、γ参数估计值也同时标记于图 5~7(b)图上方。

由图可见，负阶矩法、对数矩法2种参数估计方法得到的参数估计值几乎完全相同，拟合的NM和LM曲线几乎完全重合。而此二者又与统计得到的HS曲线几乎完全重合。这说明，这2种参数估计方法具有几乎完全相同的PDF建模精度，都能够真实反映数据的PDF特性。换言之，SαS分布S_α(γ) 可以精确反映混响数据的PDF，其参数估计选用正负矩法或对数矩法皆可。从运算速度看，对数矩法略优。下文选用对数矩法。

3.2 建模结果分析

选取(0411三亚海试)中不同方位通道、不同距离处的Ⅰ、Ⅱ、Ⅲ 3种发射信号引起的混响数据各5 430段，以对数矩法进行SαS分布建模。由于离差估值$\hat \gamma$ 仅反映数据值偏离中值或均值的程度，亦即基本反映数据功率的大小，而与数据的非高斯性关系不大，故这里仅考察特征指数估值$\hat \alpha$ 的统计规律。

表 1给出了这些数据的建模结果。3组共计16 290段数据的$ \hat \alpha$ 值皆分布于1.1~2.5，将1.1~2.0间隔0.1长度分为1段，得9个区间，由小到大表征依次削弱的非高斯性，而≥2.0则合为1个区间，表征高斯属性。表中第1列即显示了这种区间划分；第2~4列统计出各组数据$ \hat \alpha$ 值分别位于10个区间的百分比，统计的基数都是5 430段；第5列则是3组综合的百分比，统计基数为16 290段。

由表可见，各组数据的$ \hat \alpha$ 没有低于1.1的值，少部分位于1.1~1.6，多数位于1.6~2.0，各有30%左右$ \hat \alpha$ ≥2.0。考虑到时间平均量代替了统计平均量以及样本量的不足，可以把$ \hat \alpha$≥1.9的情况视为近高斯分布，则各组占比分别为：信号Ⅰ组55.4%，信号Ⅱ组56.0%，信号Ⅲ组52.2%；把$ \hat \alpha$ < 1.9的情况视为典型非高斯分布，则各组占比分别为：信号Ⅰ组44.6%，信号Ⅱ组44.0%，信号Ⅲ组47.8%；把$ \hat \alpha$ < 1.6的情况视为强非高斯分布，则各组占比分别为：信号Ⅰ组6.0%，信号Ⅱ组1.1%，信号Ⅲ组3.2%。

因此，可以说，混响数据的SαS分布特征指数近半数位于1.6~1.9，表现出一定的非高斯性。相较而言，宽带长脉冲信号引发的混响数据非高斯性略强，而窄带短脉冲信号引发的混响数据的特征指数分布范围更广(信号Ⅰ组有少量$ \hat \alpha$ 值散布于1.1~1.3，而信号Ⅱ、Ⅲ组则无)。

4 结论

1)SαS分布可对主动声呐背景干扰进行较好的概率密度拟合；

2) 在混响限制区内，主动声呐背景干扰的特征指数α近半数集中在1.6~1.9，表现出一定的非高斯性。