基于高斯混合模型的非高斯振動疲勞頻域求解方法

朱帥康，董龍雷，官 威，王 珺，李斌潮

(1.西安交通大學 航天航空學院，西安 710000；2.西安航天動力研究所，西安 710100)

機械結構的設計階段，通常使用頻域方法估計結構的疲勞壽命。當結構的應力響應服從窄帶高斯分布時，可假設雨流計數的幅值服從Rayleigh分布，再結合Miner累計損傷準則[1]及材料的S-N曲線計算出結構的疲勞壽命。如果應力幅值為寬帶信號，Rayleigh分布假設不再適用，許多學者提出了針對寬帶載荷的算法，如Wirsching-Light(WL)方法[2]，Dirlik方法[3]，Zhao-Barker(ZB)方法[4]，Tovo-Benasciutti(TB)方法[5]等，這些方法都要求載荷為高斯隨機載荷。

然而，許多隨機載荷具有較強的非高斯特性，如風載、海浪沖擊等，這些載荷產生的響應通常也呈現出非高斯性。比如，海浪對平臺的沖擊載荷峭度值大約為5，而風載的峭度值可能超過10。另外，還有一些載荷呈現多峰分布[6]，如橋梁、路面上的車輛載荷等。在不滿足單峰高斯分布載荷的條件下使用頻域疲勞計算方法，勢必會帶來較大的計算誤差，產生錯誤的疲勞壽命估計。

針對非高斯載荷，Braccesi等[7-8]提出了修正系數法，其修正系數由非高斯分布的峭度值確定。這種方法只能用于單峰非高斯分布載荷，且要求無偏斜。Winerstein[9]建立了非線性變換模型，使用Hermiter多項式將非高斯過程轉化為高斯過程，并且假設轉換前后的功率譜密度保持不變，適用于載荷的非高斯特性不強的情形。另外，Wolfsteiner等[10-11]將一個非高斯過程分解為幾個高斯過程，該方法要求非高斯分布無偏，且計算過程涉及到了應力功率譜的前8階距，實際應用中可能會產生嚴重的計算穩定性問題。

綜上所述，在多峰非高斯載荷下，目前還沒有通用的疲勞計算方法。本文針對多峰非高斯載荷，使用高斯混合模型(Gaussian mixture model,GMM)對載荷進行描述，為了滿足計算速度與精度，使用EM算法對模型參數進行求解，并結合TB方法得到多峰非高斯載荷下的結構頻域疲勞計算方法。通過一個疲勞分析示例驗證了該方法的有效性，表明該方法有更好的適用性及精度。

1 信號的非高斯特性

對于單峰的非高斯載荷，使用峭度值及偏度值可以對其進行描述[12]。對于一個高斯過程，其峭度值與偏度值分別為3和0。通常使用峭度值、偏度值與3和0的偏離程度衡量一個信號的非高斯性，且兩者可同時存在。其中，偏斜度用來衡量信號的非對稱性，如偏度值大于零則表示有更多的值分布在均值右側，如圖1所示。而峭度值描述信號在均值附近的密集程度，如峭度值大于3則表明僅有較少的值分布在均值附近，這意味著比起高斯分布，該信號有更高的“尾部”，如圖2所示。

圖1 偏斜非高斯信號Fig.1 Skew non Gaussian signal

圖2 對稱非高斯信號(峭度值>3)Fig.2 Symmetric non Gaussian signal (kurtosis>3)

另外，一些載荷信號呈現雙峰非高斯分布[13]，如一些橋梁上既有普通車輛也有重型車輛，且有許多樣本點遠離均值，呈現出明顯的非高斯性，如圖3所示。

圖3 橋梁上呈現雙峰分布的車輛載荷Fig.3 The vehicle bi-modal load on the bridge

信號的非高斯特征會對疲勞計算產生影響，在單軸應力狀態下，一個峭度值大于3的應力分布會對結構造成更大損傷。因為比起高斯信號，峭度值大于3的信號大幅值的應力循環更多。偏度值也會對疲勞計算產生影響，但是依賴于所選取的平均應力修正方法，且對計算結果的影響不如峭度值明顯。

2 高斯混合模型及EM算法

2.1 高斯混合模型

GMM可將一組數據的概率密度函數分解為多個高斯函數概率密度的線性組合

(1)

式中：xi屬于數據點X={x1,x2,…,xN}；αj為加權系數，且滿足

(2)

由于載荷數據為一維數據，G(xi;μj,Σj)為一維高斯分布函數，表達式為

(3)

將G(xi;μj,Σj)稱為GMM的第j個高斯分量，μj和Σj分別為此高斯分量的均值與方差。

理論上可以使用GMM描述任意分布的概率密度函數。程紅偉等[14]建立了二階的GMM，由于二階高斯混合模型未知參數較少，文中使用了高斯分量的前6階矩進行模型求解。當載荷呈現出多峰分布時，需要的高斯分量隨之增多，文中的求解方法不再適用，因此引入EM算法進行參數求解。EM算法本質上為一種迭代算法，對于有多個高斯分量的GMM，使用極大似然估計往往無法獲得參數的解析解，因此通過迭代來間接獲得各參數值。

2.2 EM算法

EM算法引入一組隱變量Zi={zi1,zi2,…,zik}，代表數據屬于各個高斯分量的概率

(4)

將(X,Z)稱為完整數據，其最大似然函數為

(5)

EM算法可分為E步與M步，E步即為求期望(expectation),M步為求極大(maximization)。

E步：已知觀測數據X及估計參數θi，求其隱變量的對數期望

(6)

(7)

M步：求Q(θ|θi)的極大值

(8)

可以求得，高斯分量的權系數為

(9)

高斯分量的期望向量為

(10)

高斯分量的方差為

(11)

通過上述過程，即可實現θi→θi+1的迭代。選擇終止迭代條件，可以計算隱變量的對數期望，即式(6)中的Q值，當Q值逐漸增大趨于穩定時則認為求得的參數收斂，停止迭代。

3 GMM-TB求解方法

3.1 Tovo-Benasciutti方法

Tovo和Benasciutti在2005年提出了一種計算疲勞損傷的方法，文中將疲勞損傷強度的上限與下限線性組合，從而計算出疲勞損傷，基于TB方法得到的疲勞損傷值為

(12)

式中：C和m為材料常數，由S-N曲線確定；vp為峰值密度，由式(13)確定

(13)

ξi為單邊功率譜Gxx(f)的第i階矩，表達式為

(14)

b為權重系數

(15)

參數βi用來估計譜寬度，表達式為

(16)

Γ(·)為歐拉伽馬函數，具體表達式如下

(17)

該方法在疲勞壽命計算中由許多研究證明有良好的計算精度[15-17]。

3.2 GMM-TB求解方法

對于一個非高斯過程而言，其功率譜密度曲線無法描述其非高斯特性。由GMM可以確定多個高斯分量出現的概率，因此可以根據GMM對功率譜密度進行分解。

GMM中第i個高斯分量的方差可以表示為

(18)

式中，Gi(f)為第i個高斯分量的功率譜密度。

非高斯過程X(t)的方差可以表示為

(19)

(20)

將式(18)和式(19)代入式(20)得到

(21)

假設Gi(f)沿頻率軸與GX(f)成比例關系，則Gi(f)可由式(22)確定

Gi(f)=λiGX(f)

(22)

式中，λi為比例系數，可以聯立式(18)、式(19)、式(22)得到

(23)

將式(22)、式(23)代入式(21)即可得到非高斯過程的功率譜密度的分解形式。

使用功率譜密度(power spectral density，PSD)計算疲勞損傷時，平均應力的影響可通過式(24)修正

GiT(f)=c2(σm)·Gi(f)

(24)

式中，c(σm)為平均應力修正系數，可由Goodman修正公式確定

(25)

式中：σm為平均應力；Rm為材料的拉伸極限強度。

通過以上步驟，可以確定GiT(f)，代入式(12)，通過TB計算方法得到對應損傷值，記為Di。總損傷值可由式(26)確定

(26)

從而得到疲勞損傷計算的GMM-TB計算公式

(27)

4 示 例

4.1 信號生成

首先，需要生成一雙峰分布的非高斯隨機信號，生成過程如下，首先確定生成隨機信號的功率譜密度，然后根據傅里葉逆變換，隨機信號的相位服從[0,2π]的均勻分布，信號的非高斯性可利用零均值非線性函數得到[18]，將峭度值設為5，記為x1(t)。之后將均值取為240 MPa,再生成一組隨機信號，其峭度值為10，記為x2(t)。將兩信號疊加，即可得到服從雙峰分布的非高斯載荷信號，記為xng(t)，生成的雙峰非高斯隨機載荷如圖4所示，其功率譜密度如圖5所示。生成信號的具體特性由表1列出。該信號線性坐標及對數坐標下的概率密度分布如圖6所示，為了顯示該信號的非高斯性，圖中也顯示了相應高斯分布的概率密度曲線。

圖4 雙峰非高斯信號Fig.4 Bi-modal non Gaussian signal

圖5 雙峰非高斯信號的功率譜密度曲線Fig.5 Power spectral density of bi-modal non Gaussian signal

表1 非高斯載荷信號參數Tab.1 Non-Gaussian load signal parameters

圖6 雙峰非高斯分布信號的概率密度曲線Fig.6 PDF of bi-modal non Gaussian signal

由圖6可知，比起高斯分布，該信號有著明顯更高的“尾部”，說明其大幅值載荷出現的概率更大，從而對結構造成的損傷也更大。

4.2 壽命計算及結果分析

材料的一些屬性，如S-N曲線斜率等，也會對疲勞壽命的預測精度產生較大影響。為了對方法的驗證具有一定普遍性，選取了三種材料，其S-N曲線的斜率分別取為3、4、6。材料屬性如表2所示，C為S-N曲線材料常數，Su為材料的拉伸極限強度，k為材料的S-N曲線斜率。

表2 用于疲勞計算的材料參數Tab.2 Material parameters for fatigue life calculation

按照上述EM算法流程，使用python語言對該信號進行GMM參數估計，選取一組GMM的初始值，對參數進行迭代，表3給出了當K=4時參數α的迭代結果。可以看出，隨著迭代次數的增加，參數的估計結果逐漸平穩，Q值的變化如圖7所示，也逐漸增大直至平穩。此處選取迭代50次的參數值作為GMM的估計值，在實際應用中可以增加迭代次數直至參數的估計結果趨于穩定，從而使估計值與真實值的誤差足夠小。

表3 EM算法估計GMM參數值的數值迭代結果(K=4)Tab.3 Numerical iteration results of EM algorithm for estimating of GMM (K=4)

圖7 Q值隨迭代次數增加的變化曲線Fig.7 The Q value with the increase of iteration times

由K=2開始逐漸增大K值至10，得到不同高斯分量的GMM。圖8給出了K=2，K=4，K=6，K=8，K=10時GMM的概率密度曲線，為了更好地顯示擬合效果，以對數坐標顯示。圖9給出了不同分量的GMM在大幅值處對原信號的擬合效果，結果顯示隨著分量數的增多，擬合效果逐漸變好。當K=2時，即兩個高斯分量時，GMM無法準確描述該雙峰信號，隨著K值增大，在信號兩側，載荷大幅值處，GMM與原信號的吻合程度更好。

圖8 對數坐標下不同K值GMM的概率密度曲線Fig.8 PDF of GMM with different K values in logarithmic coordinates

圖9 大幅值區域GMM的概率密度曲線Fig.9 PDF of GMM models in large value region

得到不同K值下的高斯模型后，按照第3章中的計算方法，首先由式(22)、式(23)得到載荷功率譜密度的分解形式，再根據式(24)、式(25)進行平均應力修正，得到每個高斯分量的損傷值，由式(27)得到疲勞計算壽命。由于S-N曲線斜率變化較大，為了得到合理的壽命計算值，因此在不同的K值下，對初始信號的RMS值進行了選取。表4給出了K值由2～12下GMM的計算壽命及直接使用TB方法得到的疲勞壽命，并以雨流計數法得到的疲勞壽命作為參考[19]。圖10將壽命的計算結果以圖形顯示，其中K=1代表不使用GMM，直接使用頻域法得到的壽命結果。

表4 疲勞壽命計算結果Tab.4 The results of fatigue life calculation

圖10 不同K值下的疲勞壽命計算結果Fig.10 The results of fatigue life under different K

由計算結果可以看出，對于呈雙峰分布的非高斯載荷，直接使用頻域法得到的結果與實際值相差較大，因為頻域法僅適用于單峰高斯載荷，且未進行平均應力修正；若將載荷分解為兩個高斯過程，即K=2的GMM，考慮了載荷平均應力的影響，但無法描述峭度值對材料造成的損傷；隨著K值增大，GMM中的高斯分量增多，載荷的非高斯特性被描述出來，得到的壽命結果也與真實值更加接近,材料S-N曲線斜率較小時，高斯分量取到10，其計算精度可控制在20%以內。材料S-N曲線斜率較大時，非高斯性對疲勞計算的影響更為顯著，需要較多的高斯分量計算結果才可收斂。當繼續增大K值，疲勞壽命的計算精度提高不明顯，一方面由于當K值增大時，參數收斂變得困難，GMM參數求解的誤差較大；另一方面GMM容易出現過擬合現象，K值很大時，出現了一些偏離原數據分布的高斯分量，致使計算精度變差。因此，實際進行疲勞計算時要合理選取GMM中高斯分量的個數。

5 結 論

本文提出了一種非高斯隨機載荷下對結構進行疲勞計算的頻域方法。首先引入GMM對載荷進行描述,并使用EM算法對模型參數進行求解，隨著高斯分量的增多，GMM以準確描述單峰及多峰非高斯載荷。在此基礎上對信號的功率譜密度進行分解，并結合Tovo-Benasciutti方法推導出一種多峰非高斯載荷下的頻域疲勞計算方法。為了驗證此方法，對一個雙峰分布的非高斯載荷信號進行了疲勞分析，選取三種材料，以雨流計數法作為參考，結果表明對于雙峰非高斯載荷，使用GMM可以明顯提高計算精度，且在一定范圍內，計算精度隨高斯分量的增多而提高。非高斯性在材料S-N曲線斜率較大時對計算的影響更為顯著，材料S-N曲線斜率較小時，高斯分量取至10附近，其計算精度可控制在20%以內，材料S-N曲線斜率較大時，計算所需的高斯分量較多，同時計算精度提升明顯。另外，由于文中方法以Tovo-Benasciutti頻域疲勞計算方法為基礎，因此也適用于寬帶隨機載荷的計算，對工程問題具有一定實用價值。