基于高斯混合模型的非高斯随机振动幅值概率密度函数

程红伟， 陶俊勇， 蒋 瑜， 陈 循

(1.国防科技大学 装备综合保障技术重点实验室，长沙 410073;2.国防科技大学 机电工程与自动化学院，长沙 410073)

在工程实际中，一般基于高斯假设对机械振动信号进行统计分析和研究。但是，近来越来越多的学者注意到机械振动的高斯性假设在某些情况下是不成立的，机械随机振动的非高斯特性开始受到关注。Steinwolf[1]研究了湍流边界层在机翼蒙皮所引起振动的非高斯性。Gioffré[2]，Gurley[3]和董欣[4]分别研究了风引起的振动载荷的非高斯性。Rouillard[5]深入研究了车辆振动环境的非高斯性，并提出了一种基于计算机程序的高斯分解方法。Rychlik等[6-8]对海浪引起的振动载荷的非高斯性进行了大量的统计研究。

对于平稳非高斯振动信号，概率密度函数能够全面地反映其统计特性，它决定了非高斯信号高阶矩和高阶累积量的大小。非高斯随机信号幅值概率密度函数的数学描述，主要有Edgeworth展开法、高斯变换法和最大熵法。Harremoes[9]对比分析了最大熵法和Edgeworth展开法的优缺点。Winterstein[10]分析了Edgeworth展开法的缺点，并通过对高斯概率密度函数进行非线性变换，得到基于高斯变换法的非高斯概率密度函数。在随机信号非高斯性较强的情况下，Edgeworth展开法得到的概率密度曲线会出现负值并呈现多峰态；而基于最大熵理论的方法也会使单峰的非高斯概率密度函数呈现多峰态[9]。高斯变换法的问题在于仅适用于一定的峭度范围，并且该方法计算过程复杂[10]。另外，Steinwolf[11]提出了一种基于经验信息的高斯曲线拼接法。Rouillard[5,12]提出了基于搜寻算法的高斯混合模型。高斯拼接法和Rouillard提出的高斯混合法，都能够较好地逼近非高斯随机振动信号概率密度函数。但是，这两种方法计算量大，工程中应用困难。

综上所述，针对非高斯振动信号的幅值概率密度函数，需要提出一种既能满足计算精度要求，而计算过程又相对简单的数学模型。在本研究中，我们基于高斯混合模型的数学思路，提出了基于非高斯振动信号高阶统计量的高斯混合模型。应用该数学模型对仿真非高斯振动信号和实测非高斯振动信号的幅值概率密度进行表述，通过与经验分布的对比验证了该方法的有效性。通过与其它方法的对比，进一步验证了该方法的正确性和工程应用价值。

1 非高斯振动

从理论上来说，能够全面描述一个随机过程非高斯特性的统计量和频谱函数为：高阶矩、高阶累积量和高阶谱[13-14]。但是由于随机过程高阶统计量和高阶谱的复杂性，使其在非高斯振动领域难以开展工程应用。所以，研究人员大都借助随机变量的高阶统计量来描述平稳随机过程的非高斯性，其中最常用的是标准化3阶中心矩和标准化4阶中心矩，即偏斜度γ3和峭度γ4[1,15-17]，

(1)

(2)

其中，X为非高斯随机变量；μX和σX分别为X的均值和标准差；M3和M4分别为X的3阶中心矩和4阶中心矩，可由X的概率密度函数计算得到。需要注意的是，高斯随机变量的偏斜度恒为0，峭度恒为3。

对于零均值平稳非高斯振动，通过时域样本序列可以对其偏斜度和峭度进行估计，

(3)

(4)

其中，x(t)为随机过程X(t)的样本信号；T为样本时间长度。

Rouillard[5]和Rizzi[18]分别对实际环境中的非高斯振动信号进行了大量的统计分析，并指出在实际振动环境中，很少存在严格意义上的平稳信号，所以大多随机振动的非高斯性是由短时非平稳性造成的。然而，对于短时非平稳信号，从长时间的观测序列来看，它具有平稳信号的特性。所以，在工程中对于宏观平稳，而短时非平稳的振动信号一般都基于平稳随机过程假设对其进行处理。

2 高斯混合模型

Middleton[19]在研究通信系统中多源叠加噪声信号的幅值概率分布时提出了高斯混合模型，并在通信领域得到广泛的应用。高斯混合模型的统一表达式为，

(5)

其中fNG为非高斯概率密度函数；fi(x)为第i个高斯分量的概率密度函数；αi为第i个高斯分量的权值，0 ≤αi≤ 1，∑αi=1。一般情况下，2、3阶高斯混合模型就可以给出精度足够高的非高斯概率密度函数[20]。在本研究中，我们采用2阶高斯混合模型，

fNG(x)=αf1(x)+(1-α)f2(x)

(6)

根据Middleton提出的高斯混合模型理论，需要在全面掌握各噪声源的物理特性的前提下，以Poisson分布来确定不同高斯分量的权值αi，这在机械振动信号处理中是不能实现的。

众所周知，零均值高斯概率密度函数可以由标准差完全确定。所以零均值非高斯过程的二阶高斯混合模型可以表示为：

(7)

其中σ1和σ2分别为高斯分量1和高斯分量2的标准差；α和1-α分别为高斯分量1和高斯分量2的权值。式(7)中有三个未知量：σ1、σ2和α。

可以根据非高斯振动样本信号来估计其2、4、6阶中心矩，在零均值情况下中心矩和原点矩相同，

(8)

(9)

(10)

所以有，

(11)

将式(11)代入式(9)，得方程组：

(12)

以式(8)给出的样本估计量代替真值，则有：

(13)

对于方程组(13)，通过科学计算软件求解得到未知量α，σ1和σ2，将其代入式(7)中，得到给定非高斯振动信号的幅值概率密度函数。

3 示例

为了综合验证所提出的高斯混合模型的有效性，这里给出了两个示例的详细计算过程：

(1)峭度较低的仿真非高斯振动信号；

(2)峭度较大的实测车辆振动信号。

最后，通过对4个振动信号(两个仿真信号，两个实测信号，进行误差分析，从定量的角度来验证高斯混合模型的准确性和有效性。

3.1 仿真信号

引入一个对称分布的非高斯仿真信号，如图1所示。该仿真信号的均值为0，方差为1.1976 × 103，偏斜度为0，峭度为8.1394。图1(a)为信号的时域序列；图1(b)为功率谱密度，该振动信号是宽带非高斯信号。

图1 仿真非高斯振动信号，峭度：8.1394

将图1(a)所示的样本序列代入式(8)，得，

(14)

将上式代入方程组(13)得，

(15)

将式(15)代入式(7)，得非高斯幅值概率密度函数，

(16)

图2 仿真非高斯信号幅值概率密度曲线

图2给出了基于以下四种方法得到的概率密度曲线：① 样本序列的经验分布；② 本文所提出的高斯混合模型；③ 高斯假设分布；④ 四阶Edgeworth展开法。因为时域样本序列有足够的长度，所以认为经验分布有足够高的精度，以它作为其他数学模型的参考标准。图2(a)为线性坐标下的概率密度曲线，可以清晰地显示出分布曲线中间峰值部分的差异；图2(b)为半对数坐标下的幅值概率密度曲线，可以清晰地显示出分布曲线在尾部的差异。从图2可以看出，基于高斯假设的概率密度曲线与经验分布曲线差异很大，所以工程中基于高斯假设来处理非高斯信号将引入很大的误差；基于4阶Edgeworth展开的非高斯幅值概率密度曲线出现了负值和多峰态，Edgeworth展开法只适用于峭度值很小的非高斯振动信号；基于本文提出的高斯混合法无论是在分布曲线的尖峰附近还是尾部都能精确地逼近经验分布。

3.2 实测信号

某型运输车辆载货平台实测振动信号，如图3所示。该振动信号的均值为0，方差为38.646 2，偏斜度为0，峭度为22.971 6。图3(a)为信号的时域序列；图3(b)为功率谱密度，该振动信号具有窄带非高斯特性。

图3 车辆测量非高斯振动信号，峭度为22.9716

将图3(a)中的样本序列代入式(8)得，

(17)

将式(17)代入方程组(13)得，

(18)

将式(18)代入(7)，得非高斯幅值概率密度函数，

(19)

图4给出了3.1所述的四种方法的概率密度曲线。同样采用样本经验分布作为其他数学模型的参考标准。可以看出，随着峭度的增加，基于高斯假设和Edgeworth展开的概率密度函数的误差将会增大，而本文所提出的方法仍能够较好地逼近样本信号的经验分布。

图4 车辆非高斯信号的幅值概率密度曲线

3.3 结果分析

为了进一步分析所提出方法的准确性，这里以相对均方误差来衡量各幅值概率密度曲线对经验分布曲线的偏离程度。这里定义相对均方误差为，

(20)

其中f为基于某种模型或方法得到的非高斯概率密度函数，fEM为基于样本序列得到的经验分布。

为了充分验证方法的有效性，在3.1和3.2两个示例的基础上，又分别引入了一个峭度较低的仿真信号和某型飞机实测非高斯振动信号。分别计算了本文提出的高斯混合模型概率密度函数、高斯假设下的概率密度函数和基于Edgeworth展开的概率密度函数与样本经验分布之间的相对均方误差，如表1所示。

表1 非高斯振动幅值概率密度函数相对误差r

4 结 论

基于高斯混合模型，利用高斯随机变量高阶矩之间的定量关系，结合非高斯随机振动信号的物理特性，提出了一种求解非高斯振动幅值概率密度函数的方法，即二阶高斯混合模型分解方法。该方法的数学模型简单，物理意义明确，对高峭度的非高斯振动信号仍然适用。

通过仿真和实测非高斯振动信号验证了方法的有效性和工程适用性。基于高斯混合模型的概率密度函数模型为非高斯机械振动信号的进一步研究，如疲劳分析，减振隔振等，提供了准确的统计分析工具和重要的理论支撑。

另外，二阶高斯混合模型为非高斯振动的研究提供了重要思路，它可以进一步扩展为高阶模型或应用于非高斯振动信号的频域研究，以满足更高精度和更深层次的研究需求。

参 考 文 献

[1]Steinwolf A, Rizzi S A. Non-Gaussian analysis of turbulent boundary layer fluctuating pressure on aircraft skin panels [J]. Journal of Aircraft, 2006, 43 (6): 1662-1675.

[2]Gioffré M, Gusella V, Grigoriu M. Simulation of non-Gaussian field applied to wind pressure fluctuations [J]. Probabilistic Engineering Mechanics, 2000, 5(4): 339-346.

[3]Gurley K R, Kareem A, Tognarelli M A. Simulation of a class of non-normal random processes [J]. International Journal of Non-linear Mechanics, 1996, 31 (5): 601-617.

[4]董 欣, 叶继红. 马鞍屋盖表面面积平均风压特性研究[J]. 振动与冲击, 2011, 30 (7): 21-30.

DONG Xin, YE Ji-hong. Area-averaged wind pressure on a saddle roof [J]. Journal of Vibration and Shock, 2011, 30 (7): 21-30.

[5]Rouillard V. On the non-Gaussian nature of random vehicle vibrations [C]// WCE 2007. Proceedings of the World Congress on Engineering 2007. London, U K: WCE, 2007: 1219-1224.

[6]Rychlik I, Johannesson P, Leadbetter M R. Modeling and statistical analysis of ocean-wave data using transformed Gaussian processes [J]. Marine Structures, 1997, 10(1): 13-47.

[7]Podgórski K, Rychlik I, Machado U B. Exact distributions for apparent waves in irregular seas [J]. Ocean Engineering, 2000, 27(9): 979-1016.

[8]Butler R W, Machado U B, Rychlik I. Distribution of wave crests in a non-Gaussian sea [J]. Applied Ocean Research, 2009, 31(1): 57-64.

[9]Harremo⊇s P. Maximum entropy and the Edgeworth expansion [C]// IEEE. Proceedings of the IEEE Information theory workshop. Awaji Island, Japan: IEEE, 2005: 68-71.

[10]Winterstein S R. Nonlinear vibration models for extremes and fatigue [J]. Journal of Engineering Mechanics, 1988 114 (10): 1772-1790.

[11]Steinwolf A. Approximation and simulation of probability distributions with a variable kurtosis value [J]Computational Statistics & Data Analysis, 1996, 21 (2): 163-180.

[12]Rouillard V. The synthesis of road vehicle vibrations based on the statistical distribution of segment lengths, 5th Australasian Congress on Applied Mechanics, ACAM 2007 10-12 December 2007, Brisbane, Australia

[13]Mendel J M. Tutorial on higher-order statistics (spectra) in signal processing and system theory: theoretical results and some applications [J]. Proceedings of the IEEE, 1991, 79 (3): 278-305.

[14]Swami A, Mendel J M, Nikias C L. Higher-order spectral analysis toolbox for use with Matlab, user's guide, version 2 [EB/OL].http://www.uic.edu/classes/idsc/ids594/research/BBC/hosa.pdf. 2000.12.27/2012.12.20.

[15]Benasciutti D, Tovo R. Fatigue life assessment in non-Gaussian random loadings [J]. International Journal of Fatigue, 2006, 28 (7):733-746.

[16]蒋 瑜, 陈 循, 陶俊勇,等. 指定功率谱密度、偏斜度和峭度值下的非高斯随机过程数字模拟[J]. 系统仿真学报, 2006, 18 (5): 1127-1130.

JIANG Yu, CHEN Xun, TAO Jun-yong, et al. Numerically simulating non-Gaussian random processes with specified PSD, skewness and kurtosis [J]. Journal of System Simulation, 2006, 18 (5): 1127-1130.

[17]蒋 瑜, 陶俊勇, 王得志,等. 一种新的非高斯随机振动数值模拟方法[J]. 振动与冲击, 2012, 31 (19): 169-173.

JIANG Yu, TAO Jun-yong, WANG De-zhi et al. A novel approach for numerical simulation of a non-Gaussian random vibration [J]. Journal of Vibration and Shock, 2012, 31 (19): 169-173.

[18]Rizzi S A, Przekop A, Turner T. On the response of a nonlinear structure to high kurtosis non-Gaussian random loadings [EB/OL]. http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/ 20110013658_2011014239.pdf. 2011.7.4/ 2012.12.26.

[19]Middleton D. Non-Gaussian noise models in signal processing for telecommunications: new methods and results for class A and class B noise models [J]. IEEE Transactions on Information Theory, 1999, 45(4): 1129-1149.

[20]Vastola K. Threshold detection in narrow-band non-Gaussian noise [J]. IEEE Transactions on Communications, 1984, 32 (2): 134-139.

[21]Bendat J S, Piersol A G. 随机数据分析方法 [M]. 凌福根. 北京: 国防工业出版社, 1980.