程紅偉，陶俊勇，陳 循，蔣 瑜

(1.國(guó)防科技大學(xué) 裝備綜合保障技術(shù)重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410073；2.國(guó)防科技大學(xué) 機(jī)電工程與自動(dòng)化學(xué)院，長(zhǎng)沙 410073)

工程實(shí)際中，一般基于高斯假設(shè)對(duì)隨機(jī)振動(dòng)信號(hào)進(jìn)行統(tǒng)計(jì)分析和處理。然而越來(lái)越多的研究發(fā)現(xiàn)非高斯振動(dòng)信號(hào)在實(shí)際環(huán)境中廣泛存在[1-4]，忽略振動(dòng)信號(hào)的非高斯性往往會(huì)導(dǎo)致過(guò)大的計(jì)算偏差甚至錯(cuò)誤的結(jié)果。Steinwolf等[5]研究了湍流邊界層在機(jī)翼蒙皮所引起振動(dòng)的非高斯性。葉繼紅等[6]研究了大跨屋蓋脈動(dòng)風(fēng)壓的非高斯特性。Rouillard[7-8]分析了車輛振動(dòng)環(huán)境的非高斯性，并提出了一種基于迭代算法的高斯分解方法。Rychlik等[9-11]對(duì)海浪及其引起的振動(dòng)載荷的非高斯性進(jìn)行了大量的統(tǒng)計(jì)研究。但非高斯振動(dòng)信號(hào)的概率分布函數(shù)尚需要進(jìn)一步研究。

對(duì)于一個(gè)平穩(wěn)非高斯振動(dòng)信號(hào)，概率密度函數(shù)(PDF)能夠全面反映其統(tǒng)計(jì)特性，并決定了非高斯信號(hào)的分布規(guī)律、各階統(tǒng)計(jì)矩和累積量。根據(jù)PDF曲線的對(duì)稱性可以將非高斯振動(dòng)信號(hào)分為對(duì)稱非高斯信號(hào)和偏斜非高斯信號(hào)。對(duì)稱非高斯信號(hào)的奇數(shù)階中心矩為零，偏斜非高斯信號(hào)的奇數(shù)階中心矩不為零。這兩種非高斯信號(hào)廣泛存在于實(shí)際環(huán)境中，其中對(duì)稱非高斯振動(dòng)信號(hào)多見于車輛結(jié)構(gòu)，偏斜非高斯信號(hào)多見于風(fēng)壓或波浪引起的結(jié)構(gòu)振動(dòng)。對(duì)稱非高斯信號(hào)是偏斜非高斯信號(hào)的一種特例。本文主要對(duì)典型偏斜非高斯振動(dòng)信號(hào)的概率密度函數(shù)進(jìn)行研究。非高斯隨機(jī)過(guò)程幅值概率密度函數(shù)的數(shù)學(xué)表述方法主要有Edgeworth展開法、高斯變換法和最大熵法。Harremo?s[12]對(duì)比分析了最大熵法和Edgeworth展開法的優(yōu)缺點(diǎn)。Winterstein[13]指出了Edgeworth展開法存在的問(wèn)題，并基于Hermite多項(xiàng)式展開理論提出了Winterstein變換模型。當(dāng)隨機(jī)信號(hào)非高斯性較強(qiáng)時(shí)，Edgeworth展開法得到的概率密度曲線會(huì)出現(xiàn)負(fù)值并呈現(xiàn)多峰態(tài)；采用最大熵理論的方法同樣會(huì)出現(xiàn)多峰值問(wèn)題[12]。Winterstein模型計(jì)算過(guò)程比較復(fù)雜[13]。另外，Steinwolf[14]提出了一種基于經(jīng)驗(yàn)信息的高斯曲線拼接法，該方法適用于對(duì)稱非高斯信號(hào)。

綜上所述，針對(duì)偏斜非高斯振動(dòng)信號(hào)的幅值概率密度函數(shù)，需要提出一種既能滿足計(jì)算精度要求，又相對(duì)簡(jiǎn)單的數(shù)學(xué)模型。本文基于高斯混合模型，提出了一種適用于偏斜非高斯振動(dòng)信號(hào)幅值概率密度函數(shù)的數(shù)學(xué)模型。應(yīng)用該數(shù)學(xué)模型對(duì)仿真非高斯振動(dòng)信號(hào)和實(shí)測(cè)非高斯振動(dòng)信號(hào)的幅值概率密度進(jìn)行描述，通過(guò)與經(jīng)驗(yàn)分布曲線和其他方法進(jìn)行對(duì)比分析，驗(yàn)證了該方法的有效性和較高的計(jì)算精度。

1 非高斯統(tǒng)計(jì)量

理論上講，能全面描述隨機(jī)過(guò)程非高斯特性的統(tǒng)計(jì)量為高階矩Mn(τ1,…,τn-1)或高階累積量Cn(τ1,…,τn-1)[15]。但無(wú)論是高階矩還是高階累積量，都是時(shí)間間隔變量τi的多元函數(shù)，其估計(jì)、表述和應(yīng)用都十分復(fù)雜，許多環(huán)節(jié)都可能引入較大的計(jì)算誤差。高階統(tǒng)計(jì)量的復(fù)雜性使其在非高斯振動(dòng)的定量分析中應(yīng)用較少，而主要應(yīng)用于定性分析和特征識(shí)別中。通常用靜態(tài)三階和四階統(tǒng)計(jì)量來(lái)描述平穩(wěn)隨機(jī)過(guò)程的非高斯性，即偏斜度γ3和峭度γ4[3,16]，

(1)

(2)

式中，X為非高斯隨機(jī)過(guò)程；μX和σX分別為X的均值和標(biāo)準(zhǔn)差；M3和M4分別為X的3階和4階中心矩，可由X的概率密度函數(shù)得到。高斯隨機(jī)變量的偏斜度γ3= 0，峭度γ4=3。

對(duì)于零均值平穩(wěn)非高斯振動(dòng)，通過(guò)時(shí)域樣本序列可以對(duì)其偏斜度和峭度進(jìn)行估計(jì)

(3)

(4)

式中x(t)為隨機(jī)過(guò)程X(t)的樣本序列；T為樣本的時(shí)間長(zhǎng)度。由于忽略了隨機(jī)信號(hào)的時(shí)間相關(guān)性和四階以上的高階統(tǒng)計(jì)量，偏度和峭度不能完全表示隨機(jī)振動(dòng)信號(hào)的非高斯性。但在工程中，一般考慮偏度和峭度的計(jì)算結(jié)果能夠滿足精度要求。

2 非高斯隨機(jī)振動(dòng)的高斯混合模型

2.1 高斯混合模型

Middleton[17]在研究通信系統(tǒng)多源疊加噪聲信號(hào)的幅值概率分布時(shí)提出了高斯混合模型，并在通信領(lǐng)域得到廣泛的應(yīng)用。高斯混合模型的統(tǒng)一表達(dá)式為

(5)

其中fNG為非高斯概率密度函數(shù)；fi(x)為第i個(gè)高斯分量的概率密度函數(shù)；αi為第i個(gè)高斯分量的權(quán)值，0 ≤αi≤ 1，∑αi=1。一般情況下，二階或三階高斯混合模型就可以給出精度足夠高的結(jié)果。在本研究中，我們采用二階高斯混合模型

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

(6)

2.2 高階矩分解

由于均值對(duì)概率分布的影響可以通過(guò)坐標(biāo)軸平移變換得到，因此本研究基于零均值假設(shè)展開。對(duì)于零均值偏斜非高斯過(guò)程，假設(shè)其高斯混合模型為：

(7)

式中α，m1和σ1分別為高斯分量1的權(quán)重因子、均值和標(biāo)準(zhǔn)差；1-α，m2和σ2分別為高斯分量2的權(quán)重因子、均值和標(biāo)準(zhǔn)差。均值參數(shù)m1和m2的引入使高斯混合模型可以擬合偏斜概率密度曲線。式(7)中有5個(gè)未知參數(shù)，m1，m2，σ1，σ2和α。對(duì)于零均值非高斯過(guò)程其中心矩等于原點(diǎn)矩，以下統(tǒng)稱為矩。非高斯隨機(jī)過(guò)程的一階矩為零，如式(8)

(8)

非高斯過(guò)程的二階矩為方差，如式(9)

(9)

式中Ψ1(x)和Ψ2(x)分別為高斯分量1和高斯分量2的均方值，是均值和方差的函數(shù)：

(10)

將式(10)代入式(9)，非高斯隨機(jī)過(guò)程的二階矩可以展開為

(11)

同理，非高斯過(guò)程的三階矩為

(12)

(13)

代入式(12)，非高斯隨機(jī)過(guò)程的三階矩為

(14)

類似地，非高斯隨機(jī)過(guò)程的四階矩和五階矩如式(15)和(16)所示

(15)

(16)

對(duì)于式(15)，有

(17)

將式(17)代入式(15)，則非高斯隨機(jī)過(guò)程的四階矩為

(18)

對(duì)于式(16)，有：

將式(19)代入式(16)，則非高斯隨機(jī)過(guò)程的五階矩為

(20)

實(shí)際問(wèn)題中，非高斯隨機(jī)過(guò)程的各階矩是未知的，一般根據(jù)樣本記錄得到其估計(jì)值。假設(shè)零均值非高斯過(guò)程的時(shí)間樣本序列為x(t)，則其第i階矩的估計(jì)值為

(21)

用估計(jì)值代替各階矩的理論值，聯(lián)立式(8)、(11)、(14)、(18)和(20)，得到關(guān)于未知參數(shù)α,m1,m2,σ1,σ2的五元方程組：

(22)

3 仿真算例與試驗(yàn)

為綜合驗(yàn)證所提出的高斯混合模型的有效性，分別給出了以下兩個(gè)示例：① 基于非線性變換得到的仿真加速度振動(dòng)信號(hào)；② 非高斯振動(dòng)臺(tái)懸臂梁振動(dòng)試驗(yàn)測(cè)量得到的非高斯應(yīng)力響應(yīng)信號(hào)。兩種信號(hào)均為各態(tài)歷經(jīng)隨機(jī)過(guò)程。

3.1 仿真信號(hào)

首先生成圖1所示的零均值高斯信號(hào)，信號(hào)的標(biāo)準(zhǔn)差σx= 74.85 ms-2，功率譜(PSD)如圖1(b)所示。對(duì)圖1(a)所示的高斯信號(hào)進(jìn)行如下非線性變換并去除均值，得到零均值非高斯信號(hào)z0(t)，如式(23)所示：

z(t)=x(t)+0.002x2(t)

z0(t)=z(t)-mean(z)

(23)

式中:z0(t)的時(shí)間序列及功率譜如圖2(a)和圖2(b)所示。z0(t)的標(biāo)準(zhǔn)差σz0= 76.78 ms-2，偏度γ3= 0.915 0，峭度γ4= 4.196 9。

將非高斯時(shí)間序列代z0(t)入式(21)，得到各階矩的估計(jì)值：

將上述估計(jì)結(jié)果代入式(22)，得高斯混合模型各參數(shù)的估計(jì)估計(jì)結(jié)果

圖1 仿真高斯隨機(jī)振動(dòng)信號(hào)

圖2 仿真非高斯隨機(jī)振動(dòng)信號(hào)

將上述結(jié)果代入式(7)，得到圖2(a)所示非高斯時(shí)間序列的幅值概率密度函數(shù)：

(24)

基于式(24)得到的非高斯概率密度曲線如圖3所示。圖3中同時(shí)給出了基于Edgeworth方法、Winterstein模型的概率密度曲線以及基于樣本序列的經(jīng)驗(yàn)分布曲線。圖3(a)為線性坐標(biāo)，可以清晰地顯示分布曲線中間峰值部分的差異；圖3(b)為半對(duì)數(shù)坐標(biāo)，可以清晰地顯示分布曲線在尾部的差異。通過(guò)對(duì)比分析，高斯混合模型和Winterstein模型得到的結(jié)果比較理想，Edgeworth方法計(jì)算的概率密度曲線偏離經(jīng)驗(yàn)分布比較明顯，而且曲線出現(xiàn)局部起伏，對(duì)于偏度和峭度更大的非高斯信號(hào)，這些起伏會(huì)繼續(xù)發(fā)展，導(dǎo)致出現(xiàn)多峰值甚至負(fù)值。

為進(jìn)一步分析和比較各種方法的準(zhǔn)確性，這里以相對(duì)均方誤差來(lái)衡量各幅值概率密度曲線對(duì)經(jīng)驗(yàn)分布曲線的偏離程度。相對(duì)均方誤差定義為

(25)

式中f為基于某種方法得到的非高斯概率密度函數(shù)，fEM為基于樣本序列的經(jīng)驗(yàn)分布。則圖3中，對(duì)應(yīng)于Edgeworth展開的均方誤差rEdge= 2.13%，Winterstein模型的均方誤差rWinter= 1.43%，高斯混合模型的均方誤差rGM= 1.03%。綜合圖3和均方誤差結(jié)果，可以看出高斯混合模型能夠給出圖2所示仿真非高斯信號(hào)最準(zhǔn)確的概率密度函數(shù)解析表達(dá)式。

3.2 實(shí)測(cè)信號(hào)

如圖4所示的懸臂梁結(jié)構(gòu)，其材料為鋁合金2024-T3。懸臂梁的幾何尺寸見圖4，懸臂梁根部通過(guò)夾具固定在振動(dòng)臺(tái)上構(gòu)成如圖所示的基礎(chǔ)激勵(lì)振動(dòng)系統(tǒng)。基礎(chǔ)激勵(lì)的輸入信號(hào)是加速信號(hào)，如圖5所示。輸入加速度為對(duì)稱非高斯信號(hào)，偏度γ3= 0，峭度γ4= 6，標(biāo)準(zhǔn)差σ= 10 g(g表示重力加速度)。振動(dòng)試驗(yàn)是在豎直方向上進(jìn)行的，由于重力和其他非線性因素的影響，圖4所示的懸臂梁根部的應(yīng)力響應(yīng)為偏斜非高斯過(guò)程，去除均值以后的應(yīng)力信號(hào)如圖6所示，信號(hào)標(biāo)準(zhǔn)差σ=41 MPa，偏度γ3=-0.7102，峭度γ4=4.7977。將圖6所示的非高斯應(yīng)力響應(yīng)序列代入式(21)得到各階矩的估計(jì)值：

將上述結(jié)果代入式(22)，得到高斯混合模型各個(gè)參數(shù)的估計(jì)結(jié)果

圖3 仿真非高斯信號(hào)幅值概率密度曲線

圖4 懸臂梁結(jié)構(gòu)基礎(chǔ)激勵(lì)振動(dòng)試驗(yàn)(單位：mm)

圖5 非高斯振動(dòng)試驗(yàn)輸入信號(hào)

圖6 非高斯振動(dòng)試驗(yàn)的應(yīng)力響應(yīng)

將上述結(jié)果代入式(7)，得到圖6(a)所示非高斯振動(dòng)信號(hào)的幅值概率密度函數(shù)：

(26)

圖7中分別給出了基于高斯混合模型、Edgeworth展開法和Winterstein模型的概率密度曲線以及基于樣本序列的經(jīng)驗(yàn)分布曲線。根據(jù)式(26)Edgeworth展開法相對(duì)于經(jīng)驗(yàn)分布的均方誤差為rEdge=1.39%；Winterstein模型的均方誤差為rWinter=0.91%，高斯混合模型的均方誤差為rGM= 0.30%。通過(guò)圖7的分析和均方誤差的定量比較發(fā)現(xiàn)，對(duì)于圖6所示的實(shí)測(cè)非高斯應(yīng)力信號(hào)，基于高斯混合模型能夠給出其概率密度曲線的最優(yōu)解析表達(dá)式。

對(duì)于工程中的偏斜非高斯振動(dòng)信號(hào)，其偏度，-1.2<γ3<1.2。示例3.1中仿真信號(hào)的偏度為0.915 0，示例3.2中實(shí)測(cè)信號(hào)的偏度為-0.710 2。二者分別接近常見非高斯信號(hào)偏度的上限和下限?？偨Y(jié)示例的分析結(jié)果，可以得出結(jié)論：基于高斯混合模型的數(shù)學(xué)表達(dá)式能夠準(zhǔn)確地表示偏斜非高斯信號(hào)的幅值概率密度。

圖7 實(shí)測(cè)非高斯振動(dòng)應(yīng)力信號(hào)幅值概率密度曲線

4 結(jié) 論

基于高斯混合模型，利用高斯隨機(jī)變量高階中心矩和原點(diǎn)矩之間的關(guān)系，提出了一種求解偏斜非高斯振動(dòng)幅值概率密度函數(shù)的方法。該方法的數(shù)學(xué)模型簡(jiǎn)單，能夠準(zhǔn)確表示偏斜非高斯振動(dòng)信的概率密度。

通過(guò)仿真和實(shí)測(cè)非高斯振動(dòng)信號(hào)驗(yàn)證了所提出方法的有效性和工程適用性。

本文提出的基于高斯混合模型的概率密度函數(shù)模型為偏斜非高斯振動(dòng)信號(hào)的相關(guān)研究研究(如疲勞分析，減振隔振等)提供了準(zhǔn)確的統(tǒng)計(jì)分析工具和重要的理論支撐。

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

[2] 王旭，黃鵬，顧明.海邊坡角可調(diào)試驗(yàn)房風(fēng)載荷現(xiàn)場(chǎng)實(shí)測(cè)研究[J].振動(dòng)與沖擊，2012，31(5): 176-182.

WANG Xu，HUANG Peng，GU Ming.Field investigation on wind loads of a low building with adjustable roof pitch near see [J].Journal of Vibration and Shock，2012，31(5): 176-182.

[3] 何軍.非高斯載荷作用下結(jié)構(gòu)首次失效時(shí)間分析的Monte-Carlo模擬方法[J].振動(dòng)與沖擊，2007，26(3): 59-60.

HE Jun.Monte-Carlo simulation for first failure time of structures excited by non-Gaussian load [J].Journal of Vibration and Shock，2007，26(3): 59-60.

[4] Aberg S，Podgórski K，Rycklik I.Fatigue damage assessment for a spectral model of non-Gaussian random loads [J].Probabilistic Engineering Mechanics，2009，24(4): 608-617.

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

[6] 葉繼紅，侯信真.大跨屋蓋脈動(dòng)風(fēng)壓得非高斯特性研究[J].振動(dòng)與沖擊，2010，29(7): 9-15.

YE Ji-hong，HOU Xin-zhen.Non-Gaussian features of fluctuating wind pressures on long span roofs [J].Journal of Vibration and Shock，2010，29(7): 9-15.

[7] Rouillard V.The synthesis of road vehicle vibrations based on the statistical distribution of segment lengths [C]//ACAM.Proceedings of 5th Australasian Congress on Applied Mechanics.Brisbane，Australia: ACAM，2007:1-6.

[8] Rouillard V.On the non-Gaussian nature of random vehicle vibrations [C]// WCE 2007.Proceedings of the World Congress on Engineering 2007.London，UK: WCE，2007.1219-1224.

[9] 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.

[10] Podgórski K，Rychlik I，Machado U B.Exact distributions for apparent waves in irregular seas [J].Ocean Engineering，2000，27(9): 979-1016.[11] 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.

[12] 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.

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

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

[15] 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.

[16] 蔣瑜，陶俊勇，王得志,等.一種新的非高斯隨機(jī)振動(dòng)數(shù)值模擬方法[J].振動(dòng)與沖擊，2012，31(9): 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(9): 169-173.

[17] 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.