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

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

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

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

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

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

1 信號的非高斯特性

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

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

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

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

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

信號的非高斯特征會對疲勞計(jì)算產(chǎn)生影響，在單軸應(yīng)力狀態(tài)下，一個峭度值大于3的應(yīng)力分布會對結(jié)構(gòu)造成更大損傷。因?yàn)楸绕鸶咚剐盘?，峭度值大?的信號大幅值的應(yīng)力循環(huán)更多。偏度值也會對疲勞計(jì)算產(chǎn)生影響，但是依賴于所選取的平均應(yīng)力修正方法，且對計(jì)算結(jié)果的影響不如峭度值明顯。

2 高斯混合模型及EM算法

2.1 高斯混合模型

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

(1)

式中：xi屬于數(shù)據(jù)點(diǎn)X={x1,x2,…,xN}；αj為加權(quán)系數(shù)，且滿足

(2)

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

(3)

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

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

2.2 EM算法

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

(4)

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

(5)

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

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

(6)

(7)

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

(8)

可以求得，高斯分量的權(quán)系數(shù)為

(9)

高斯分量的期望向量為

(10)

高斯分量的方差為

(11)

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

3 GMM-TB求解方法

3.1 Tovo-Benasciutti方法

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

(12)

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

(13)

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

(14)

b為權(quán)重系數(shù)

(15)

參數(shù)βi用來估計(jì)譜寬度，表達(dá)式為

(16)

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

(17)

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

3.2 GMM-TB求解方法

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

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

(18)

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

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

(19)

(20)

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

(21)

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

Gi(f)=λiGX(f)

(22)

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

(23)

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

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

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

(24)

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

(25)

式中：σm為平均應(yīng)力；Rm為材料的拉伸極限強(qiáng)度。

通過以上步驟，可以確定GiT(f)，代入式(12)，通過TB計(jì)算方法得到對應(yīng)損傷值，記為Di?？倱p傷值可由式(26)確定

(26)

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

(27)

4 示 例

4.1 信號生成

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

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

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

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

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

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

4.2 壽命計(jì)算及結(jié)果分析

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

表2 用于疲勞計(jì)算的材料參數(shù)Tab.2 Material parameters for fatigue life calculation

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

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

圖7 Q值隨迭代次數(shù)增加的變化曲線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的概率密度曲線，為了更好地顯示擬合效果，以對數(shù)坐標(biāo)顯示。圖9給出了不同分量的GMM在大幅值處對原信號的擬合效果，結(jié)果顯示隨著分量數(shù)的增多，擬合效果逐漸變好。當(dāng)K=2時，即兩個高斯分量時，GMM無法準(zhǔn)確描述該雙峰信號，隨著K值增大，在信號兩側(cè)，載荷大幅值處，GMM與原信號的吻合程度更好。

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

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

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

表4 疲勞壽命計(jì)算結(jié)果Tab.4 The results of fatigue life calculation

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

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

5 結(jié) 論

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