0 引　言

在主动声呐信号处理领域中，海洋混响是最常见的背景干扰。传统信号处理方法通常将混响视为局部平稳的高斯色噪声，然后对混响进行预白化处理，以期将混响干扰转变为高斯白噪声干扰^{[1 - 2]}。但随着声呐向低频、大孔径和高分辨发展，混响的非高斯特性逐渐凸显，而对混响进行非高斯概率密度建模已成为一种趋势^{[3 − 7]}。

混合高斯（Gaussian Mixture，GM）分布模型^[5]和对称α稳定分布（Symmetric alpha-stable Distribution，SαS）模型 ^{[6 - 7]}是比较常用的2种非高斯模型，下文简称GM分布模型和SαS分布模型。在对海洋混响进行概率密度建模时，往往是将混响数据按时间先后顺序分成若干段，再对每一段数据预白、建模。但由于海洋混响局部平稳，每一段混响数据的概率密度特性都有所不同，这会导致某些段数据采用GM分布模型建模效果较好，但某些段数据采用SαS分布模型建模更为适合，即很难用一个确定的分布模型始终准确对混响建模，因此无论是GM分布模型还是SαS分布模型，其建模的稳健性都较弱，都会发生模型失配的情况。一旦模型失配，可能会导致后续信号处理效果明显下降^[8]。因此，需要寻找一种稳健的概率密度建模方法。

本文参考GM分布模型和SαS分布模型的特点，采用高斯-污染混合分布模型对海洋混响建模，基于最大熵原则确定污染分布类型，并在此基础上进行参数估计，最后基于实际混响数据验证了所提模型的稳健性。

1 混响建模常用非高斯模型

对海洋混响进行非高斯建模时，GM分布模型和SαS分布模型是2种较为常见的概率分布模型。

1.1 GM分布模型

GM分布模型的概率密度函数（Probability Density Function，PDF）一般可表示为：

$ {f_{gm}}\left( x \right) = \sum\limits_{m = 1}^M {{\varepsilon _m}} {\varphi _m}\left( x \right)。$

(1)

式中：$ M $为混合高斯模型的阶数；$ {\varphi _m}\left( x \right) $为高斯分布$ N\left( {{\mu _m},\sigma _m^2} \right) $的PDF；$ {\varepsilon _m}\left( {m = 1,2,...,M} \right) $称为混合参数，且满足$ {\varepsilon _m} > 0,\sum\limits_{m = 1}^M {{\varepsilon _m}} = 1 $。

对混响建模时，多采用二阶零均值混合高斯模型，其PDF为：

$ {f_{gm}}\left( x \right) = \left( {1{\text{ - }}\varepsilon } \right){\varphi _1}\left( x \right) + \varepsilon {\varphi _2}\left( x \right)。$

(2)

式中：$ {\varphi _1}\left( x \right) = \displaystyle\frac{1}{{\sqrt {2{\text{π }}} {\sigma _1}}}{{\text{e}}^{ - \frac{{{x^2}}}{{2\sigma _1^2}}}} $、$ {\varphi _2}\left( x \right) = \displaystyle\frac{1}{{\sqrt {2{\text{π }}} {\sigma _2}}}{{\text{e}}^{ - \frac{{{x^2}}}{{2\sigma _2^2}}}} $分布为高斯分布$ N\left( {0,\sigma _1^2} \right) $、$ N\left( {0,\sigma _2^2} \right) $的PDF。

1.2 SαS分布模型

SαS分布模型没有显式的PDF表达式，一般是通过对其特征函数进行傅里叶变换得到其PDF，其表达式如下：

$ {f_\alpha }\left( x \right) = \frac{1}{{\text{π}} }\int_0^\infty {\exp \left( { - \gamma {t^\alpha }} \right)\cos \left( {xt} \right){{\rm{d}}} t}。$

(3)

式中：$\alpha \in (0,2]$为特征指数，其值越小，表示混响的拖尾越厚，非高斯特性越显著；$\gamma \in (0,\infty )$为尺度参数，表征混响偏离其均值的程度。

式(3)一般没有封闭的解析表达式，需要以数值积分方式计算，这对其应用造成很大不便。LI等 ^[9]指出，在$1 \leqslant \alpha \leqslant 2$时，该分布的PDF可由式(4)近似表示。

$ {f_\alpha }\left( x \right) = \frac{{1 - \varepsilon }}{{2{\sigma _\alpha }\sqrt {\text {π}} }}\exp \left( { - \frac{{{x^2}}}{{4\sigma _\alpha ^2}}} \right) + \frac{{\varepsilon {\sigma _\alpha }}}{{{\text {π}} \left[ {{x^2} + \sigma _\alpha ^2} \right]}}。$

(4)

式中：${\sigma _\alpha } = {\gamma ^{1/\alpha }}$称为离差；$\varepsilon \left( {0 \leqslant \varepsilon \leqslant 1} \right)$称为混合参数，它仅与$\alpha $值有关，且$\varepsilon \approx \displaystyle\frac{{4 - {\alpha ^2}}}{{3{\alpha ^2}}}$。

事实上，从式(4)可看出，它是均值为0、方差为$4\sigma _\alpha ^2$的高斯分布和尺度参数为${\sigma _\alpha }$的柯西分布混合而成，也即，SαS分布模型可由高斯分布和柯西分布的组合进行近似。

2 高斯-污染混合分布模型 2.1 高斯-污染混合分布模型

高斯-污染混合分布模型的PDF表达式如下^[10]：

$ f\left( x \right) = \left( {1 - \varepsilon } \right)\varphi \left( x \right) + \varepsilon g\left( x \right)。$

(5)

式中：$ \varphi \left( x \right) $为高斯分布$ N\left( {\mu ,{\sigma ^2}} \right) $的PDF；$ g\left( x \right) $为污染分布的PDF；$\varepsilon $称为污染系数。

从式(5)可看出，高斯-污染混合分布是高斯分布和一个污染分布的线性组合。当污染分布已知时，它是一个确知的非高斯分布，当污染分布类型未知时，它是一个未知分布，但部分统计特性已知。

另外，对比式(2)、式(4)和式(5)可知，通过取不同的$ \varphi \left( x \right) $和$ g\left( x \right) $，二阶零均值混合高斯模型和SαS分布模型（其近似形式）均可视为高斯-污染分布模型的特例。

由于混响是零均值的，不妨假设高斯分布和污染分布的均值都为0，即$ \varphi \left( x \right) $为高斯分布$ N\left( {0,{\sigma ^2}} \right) $的PDF，此时式(5)可进一步写为：

$ f\left( x \right) = \frac{{1 - \varepsilon }}{{\sqrt {2{\text{π}} {\sigma ^2}} }}\exp \left( { - \frac{{{x^2}}}{{2{\sigma ^2}}}} \right) + \varepsilon g(x)。$

(6)

式中：$ g\left( x \right) $所表征的污染分布模型其均值为0，即满足$ \int_{ - \infty }^{ + \infty } {xg\left( x \right)} = 0 $。

2.2 最大熵污染分布

当污染分布未知时，高斯-污染混合分布模型是一个未知的分布模型，因此对混响建模时首先需要确定污染分布模型。但满足均值为0的概率分布模型有无穷多个，需要选择一个稳健的污染分布类型，以适应混响的时变特性。所谓稳健，是指无论混响的真实概率密度特性如何变换，建模模型都不会与之产生较大偏差。这与最大熵原则比较接近，可基于最大熵原则确定污染分布模型。

最大熵原则^[11]是由Jaynes在1957年提出，即：如果一种概率分布能在附加约束条件下使得其信息熵最大，则它是最少偏见的，也即最可能出现的分布。

对于污染分布模型$ g\left( x \right) $，其信息熵^[11]定义为：

$ H = - \int_R {g\left( x \right)\ln g\left( x \right){\rm{d}}x}。$

(7)

式中：$ R $为变量$ x $的取值边界。

在给定约束条件下使得信息熵最大的污染分布模型即为最大熵污染分布模型。通常选择整数阶样本原点矩作为约束条件，由此构建的求解最大熵污染分布模型的优化问题为：

$ \begin{gathered} \max H = - \int_R {g\left( x \right)\ln g\left( x \right){\rm{d}}x}，\\ {\rm{s}}.{\rm{t}}.\;\;\;\;\int_R {g\left( x \right){\rm{d}}x = 1}，\\ \;\;\;\;\;\;\;\;\int_R {{x^i}g\left( x \right){\rm{d}}x} = {\mu _i},\;i = 1,2,...,N。\\ \end{gathered} $

(8)

式中：$ {\mu _i} $为第$ i $阶样本原点矩；$ N $为采用的最高的样本原点矩的阶数；$ \mathrm{s}.\mathrm{t}. $为约束符号。

样本原点矩作为总体矩的替代，包含了样本的统计信息。一阶矩反映样本的均值信息；三阶矩反映样本的偏态信息。理论上，约束的矩阶数越多，概率密度建模越准确，但对海洋混响进行概率密度建模时，选择的概率密度模型为高斯-污染混合分布，因此样本的二阶、三阶矩反映的是高斯-污染混合分布的二阶矩和三阶矩，并不是污染分布的二阶矩和三阶矩。对于污染分布而言，已知的信息只有其均值为0，也即一阶矩为0，因此式(8)所示的优化问题可简化为：

$ \begin{gathered} \max H = - \int_R {g\left( x \right)\ln g\left( x \right){\rm{d}}x} ，\\ {\rm{s}}.{\rm{t}}.\;\;\;\;\int_R {g\left( x \right){\rm{d}}x = 1} ，\\ \;\;\;\;\;\;\;\;\int_R {xg\left( x \right){\rm{d}}x} = 0。\\ \end{gathered} $

(9)

对于式(9)所示的约束优化问题，可采用拉格朗日乘子式将其转化为无约束优化问题，在信息熵最大时，$ g\left( x \right) $的解析式为^[12]：

$ g\left( x \right) = \exp \left( { - {\lambda _0} - {\lambda _1}\left| x \right|} \right)。$

(10)

式中：$ {\lambda }_{0}、{\lambda }_{1} $为模型参数。将式(9)中的约束条件代入式(10)中，可得$ {e^{{\lambda _0}}} = 2/{\lambda _1} $，并令$ b = 1/{\lambda _1} $，代入式(10)可得：

$ g\left( x \right) = \frac{1}{{2b}}{e^{ - \frac{{\left| x \right|}}{b}}} 。$

(11)

式(11)为拉普拉斯分布的概率密度函数，$ b $又称为拉普拉斯分布的形状参数。

由上述推导可知，均值约束条件下最大熵高斯-污染分布混合模型即为高斯-拉普拉斯混合分布模型（Gaussian-Laplace mixture, GL），下文简称GL分布模型。将式(11)代入式(6)，可得GL分布模型的PDF为：

$ f\left( x \right) = \frac{{1 - \varepsilon }}{{\sqrt {2\pi {\sigma ^2}} }}\exp \left( { - \frac{{{x^2}}}{{2{\sigma ^2}}}} \right) + \frac{\varepsilon }{{2b}}{e^{ - \frac{{\left| x \right|}}{b}}} 。$

(12)

3 GL分布模型参数估计

记GL分布模型的概率密度函数为${f_{gl}}\left( x \right)$，将式(12)重写为：

$ {f_{gl}}\left( x \right) = \left( {1 - \varepsilon } \right){f_g}\left( x \right) + \varepsilon {f_l}\left( x \right)。$

(13)

式中：$ {f_g}\left( x \right) = \displaystyle\frac{1}{{\sqrt {2{\text{π }}} \sigma }}{{{e}}^{ - \frac{{{x^2}}}{{2{\sigma ^2}}}}} $为混合分布中高斯分量的PDF；$ {f_l}\left( x \right) = \displaystyle\frac{1}{{2b}}{e^{ - \frac{{\left| x \right|}}{b}}} $为混合分布中拉普拉斯分量的PDF。

由式(13)可知，GL分布模型由高斯分布方差$ {\sigma ^2} $、拉普拉斯分布形状参数$ b $以及污染系数$ \varepsilon $所组成的参数组${\mathbf{\theta }} = {[\varepsilon ,{\sigma ^2},b]^{\text{T}}}$唯一确定。对混响进行GL分布模型建模的过程，也即根据样本序列$ {{X}} = \{ {x_1},{x_2},...,{x_{\text{N}}}\} $，估计PDF参数组${\mathbf{\theta }} = {[\varepsilon ,{\sigma ^2},b]^{\text{T}}}$的过程。

3.1 EM算法估计步骤

EM（Expectation Maximization）算法^[13]，又称期望最大化算法，在GM分布模型参数估计中应用广泛。由于GL分布模型也属于混合分布模型的一种，故对EM算法进行改进，将其应用于GL分布模型的参数估计中。具体推导过程不再详述，直接给出EM算法估计GL分布模型参数的步骤。

步骤1　根据式(14)给定GL分布模型PDF参数初值${{\mathbf{\theta }}_0} = {[{\varepsilon _0},\sigma _0^2,{b_0}]^{\text{T}}}$

$ \left\{ \begin{gathered} {\varepsilon _0} = 0.5，\\ \sigma _0^2 = \sqrt {\frac{{Var}}{2}} ，\\ {b_0} = \frac{{\sqrt {Var} }}{2}。\\ \end{gathered} \right. $

(14)

式中：$ Var $为混响数据方差。式(14)本质上是按高斯分量和拉普拉斯分量各占一半确定的参数初值。如果初值设置不当，可能导致估计结果错误。

步骤2　根据式(15)计算$ {x_i}\left( {i = 1,2,...N} \right) $属于第$ k $个分量的后验概率$ {{P}}(k\left| {{x_i}} \right.) $

$\left\{ \begin{aligned} &{{P}}\left( {k = 1|{x_i}} \right) = \frac{{\left( {1 - \varepsilon } \right){f_g}\left( {{x_i}} \right)}}{{\left( {1 - \varepsilon } \right){f_g}\left( x \right) + \varepsilon {f_l}\left( x \right)}}，\\ &{{P}}\left( {k = 2|{x_i}} \right) = \frac{{\varepsilon {f_l}\left( x \right)}}{{\left( {1 - \varepsilon } \right){f_g}\left( x \right) + \varepsilon {f_l}\left( x \right)}}。\\ \end{aligned}\right. $

(15)

式中：$ k = 1 $为高斯分量；$ k = 2 $为拉普拉斯分量。

步骤3　根据式(16)对参数组${\mathbf{\theta }} = {[\varepsilon ,{\sigma ^2},b]^{\text{T}}}$进行迭代更新，得到$\hat {\mathbf{\theta }} = {[ \varepsilon ,{{\sigma ^2}},\hat b]^{\text{T}}}$。

$ \left\{ \begin{gathered} \hat \varepsilon = \frac{{\displaystyle\sum\limits_{i = 1}^N {{{P}}\left( {k = 2|{x_i}} \right)} }}{{\displaystyle\sum\limits_{i = 1}^N {{{P}}\left( {k = 1|{x_i}} \right) + \displaystyle\sum\limits_{i = 1}^N {{\text{P}}\left( {k = 2|{x_i}} \right)} } }}，\\ {{\sigma ^2}} = \frac{{\displaystyle\sum\limits_{i = 1}^N {\left( {{{P}}\left( {k = 1|{x_i}} \right)x_i^2} \right)} }}{{\displaystyle\sum\limits_{i = 1}^N {{{P}}\left( {k = 1|{x_i}} \right)} }} ，\\ \hat b = \frac{{\displaystyle\sum\limits_{i = 1}^N {\left( {{{P}}\left( {k = 2|{x_i}} \right){x_i}} \right)} }}{{\displaystyle\sum\limits_{i = 1}^N {{{P}}\left( {k = 2|{x_i}} \right)} }} 。\\ \end{gathered} \right. $

(16)

步骤4　当前后2次估计值的差值平方和小于给定值时，迭代终止。

EM算法估计GL分布模型参数的具体流程图如图1所示。

3.2 仿真实例

取GL分布模型的PDF参数组为θ=[ε=0.3, σ²=1, b=2]^T，生成10000个样本序列，如图2所示。然后按照图1所示的EM估计算法流程进行参数估计，估计结果为$\hat {\mathbf{\theta }} = {[\varepsilon {\text{ = }}0.30, {{\sigma ^2}}{\text{ = }}0.99,\hat b{\text{ = }}1.95]^{\text{T}}}$，估计值与理论值非常接近。

样本的真实概率密度函数曲线和估计得到的概率密度函数曲线对比如图3所示，可以明显看出估计得到的曲线和理论值基本一致，这表明EM算法用以估计GL分布模型的PDF参数可行。

4 海洋混响建模分析

某次海试时（实验海区深度为300 m，泥沙底质，拖曳声源和拖曳线列阵平行布放，深度相同），发射的主动脉冲为LFM信号，脉冲宽度为0.2 s，中心频率为1500 Hz，带宽为500 Hz，采样频率为5000 Hz。图4(a)为该脉冲所激励的一段混响信号，图4(b)为混响信号预白后结果。

对预白后的混响进行PDF建模，建模时分段进行，每段数据长度为0.2 s，步长为0.1 s（相当于重合率50%，总段数为40）。对每一段混响数据分别按照高斯分布模型、GM分布模型、SαS分布模型以及GL分布模型建模。由于混响的真实概率分布未知，故采用经典的核密度估计结果作为混响的真实概率分布。

图5分别给出了第8段数据和第10段数据的建模结果。整体来看，高斯分布模型建模结果与理论值偏差较大，表明混响数据呈现出非高斯特性。其次，从图5(a)可知，对于第8段混响数据，GL分布模型和SαS分布模型建模结果与理论值偏差较小；而对于第10段混响数据，从图5(b)可知GL分布模型和GM分布模型建模结果与理论值偏差较小。图5表明，对于不同的数据段，混响的概率密度特性确实有所不同，并且SαS分布模型和GM分布模型各有适用，但GL分布模型建模结果与真实分布偏差始终较小，表现出建模的稳健性。

为进一步说明GL分布模型建模的稳健性，采用Kullback-Leibler (KL) 散度^[14]定量描述建模模型与混响真实概率密度模型的偏差。KL散度的表达式如下：

$ {D_{KL}}(p||q) = \int_{ - \infty }^{ + \infty } {p(x)\ln \frac{{p(x)}}{{q(x)}}{\rm{d}}x} 。$

(17)

式中：$p\left( x \right)$为混响的真实概率密度函数；$q\left( x \right)$为建模得到的概率密度函数。KL散度值越高，表示建模模型与混响真实概率密度之间的偏差越大，反之，偏差越小。

不同建模模型的KL散度对比如图6所示。需要说明的是，由于混响的真实概率密度模型未知，故采用经典的核密度估计结果作为真实的概率密度值。

从图6中首先可知，采用高斯分布模型建模的KL散度值普遍较高，说明混响的非高斯性较强。其次，采用GM分布模型和SαS分布模型建模，其整体KL散度值相比高斯分布模型建模低，但起伏较大，有的数据段GM分布模型的KL散度值低，有的数据段则SαS分布模型的KL散度值低，这表明无论GM分布模型还是SαS分布模型，都无法做到始终准确对海洋混响建模，稳健性较弱。最后可知，采用GL分布模型建模，不仅KL散度曲线较为平滑，且取值始终较低，这表明GL分布模型建模与混响真实概率密度的偏差始终较小，稳健性强。

5 结　语

本文采用高斯-污染混合分布模型对混响建模，即将混响的概率密度模型视为高斯分布和未知的污染分布的线性组合，然后基于最大熵原则确定了污染分布模型为拉普拉斯分布模型，然后应用EM算法估计GL分布模型的概率密度参数。最后基于真实海洋混响数据建模仿真，表明应用GL分布模型对海洋混响进行概率密度建模具有稳健性，对后续的信号检测、参数估计等处理具有重要意义。