偏斜非高斯隨機(jī)振動信號幅值概率密度函數(shù)研究

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

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

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

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

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

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

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

(1)

(2)

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

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

(3)

(4)

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

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

2.1 高斯混合模型

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

(5)

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

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

(6)

2.2 高階矩分解

由于均值對概率分布的影響可以通過坐標(biāo)軸平移變換得到，因此本研究基于零均值假設(shè)展開。對于零均值偏斜非高斯過程，假設(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個未知參數(shù)，m1，m2，σ1，σ2和α。對于零均值非高斯過程其中心矩等于原點(diǎn)矩，以下統(tǒng)稱為矩。非高斯隨機(jī)過程的一階矩為零，如式(8)

(8)

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

(9)

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

(10)

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

(11)

同理，非高斯過程的三階矩為

(12)

(13)

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

(14)

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

(15)

(16)

對于式(15)，有

(17)

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

(18)

對于式(16)，有：

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

(20)

實(shí)際問題中，非高斯隨機(jī)過程的各階矩是未知的，一般根據(jù)樣本記錄得到其估計(jì)值。假設(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)證所提出的高斯混合模型的有效性，分別給出了以下兩個示例：① 基于非線性變換得到的仿真加速度振動信號；② 非高斯振動臺懸臂梁振動試驗(yàn)測量得到的非高斯應(yīng)力響應(yīng)信號。兩種信號均為各態(tài)歷經(jīng)隨機(jī)過程。

3.1 仿真信號

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

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

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

(23)

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

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

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

圖1 仿真高斯隨機(jī)振動信號

圖2 仿真非高斯隨機(jī)振動信號

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

(24)

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

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

(25)

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

3.2 實(shí)測信號

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

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

圖3 仿真非高斯信號幅值概率密度曲線

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

圖5 非高斯振動試驗(yàn)輸入信號

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

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

(26)

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

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

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

4 結(jié) 論

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

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

本文提出的基于高斯混合模型的概率密度函數(shù)模型為偏斜非高斯振動信號的相關(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)場實(shí)測研究[J].振動與沖擊，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)首次失效時間分析的Monte-Carlo模擬方法[J].振動與沖擊，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] 葉繼紅，侯信真.大跨屋蓋脈動風(fēng)壓得非高斯特性研究[J].振動與沖擊，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ī)振動數(shù)值模擬方法[J].振動與沖擊，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.