鄭國峰，陳柏先，隗寒冰，楊昊民，嚴璐瑤

（1.重慶交通大學機電與車輛工程學院，重慶 400074；2.中國汽車工程研究院股份有限公司，重慶 401122）

0 引言

全壽命周期內(nèi)的服役載荷譜，對零部件疲勞耐久性的準確開發(fā)具有重要意義。但對零部件全壽命周期下部件的服役載荷進行采集，需要耗費大量成本?？紤]到零部件在典型工況下的服役載荷具有相應(yīng)的統(tǒng)計規(guī)律，可通過掌握部分隨機載荷統(tǒng)計規(guī)律，再基于統(tǒng)計規(guī)律進行外推，實現(xiàn)對未采集載荷譜的掌握。載荷譜外推的本質(zhì)就是通過采集載荷譜的統(tǒng)計規(guī)律預(yù)測相同工況下未采集載荷譜。

時域外推法認為載荷譜時間序列中超過閾值的極大值和極小值服從某種分布，基于概率密度分布函數(shù)可以實現(xiàn)對時間序列規(guī)律的掌握和預(yù)測。楊子涵等［1］提出基于經(jīng)驗?zāi)B(tài)分解－峰值超閾值模型的載荷譜外推方法，通過對外推負荷篩選、閾值限定和極值排序，解決了傳統(tǒng)時域外推方法在載荷外推時結(jié)果失真的問題。He等［2］對數(shù)控機床時域載荷外推方法中閾值的求解問題，提出了一種基于灰色關(guān)聯(lián)分析的最優(yōu)閾值確定模型，進一步提高數(shù)控機床載荷譜精度。鄭國峰等［3］對廣義Pareto分布模型形狀參數(shù)的2種情況進行研究，發(fā)現(xiàn)形狀參數(shù)為零時的載荷外推結(jié)果比形狀參數(shù)不為零時的結(jié)果更為保守。Wang等［4］針對軌道車輛動態(tài)應(yīng)力數(shù)據(jù)進行研究，得出一種軌道車輛長時間服役條件下疲勞壽命的評估方法，即通過極值理論對動態(tài)應(yīng)力進行外推，發(fā)現(xiàn)當采集載荷過小時，會造成一定的誤差。陳晉市等［5］提出將主成分分析與多種擬合檢驗準則相結(jié)合選取最優(yōu)閾值的載荷譜外推方法，用于挖掘機主泵載荷外推研究，但受主成分分析方法影響較大。Yang等［6］提出一種改進的閾值選擇方法，即分析極值樣本的平均值來獲取最佳閾值對拖拉機時域載荷進行外推，其結(jié)果具有較高的精度，同時保留了載荷譜的時頻特性?；谝陨蠁栴}，本文以均方誤差為閾值選擇標準，對商用車駕駛室穩(wěn)定桿扭轉(zhuǎn)載荷進行時域外推，該方法受采集載荷樣本的影響較小，同時能夠保留載荷譜的統(tǒng)計特性。對于時域外推法中載荷譜的閾值選取，常用的選取方法有經(jīng)驗法、灰色關(guān)聯(lián)度法［2，7－8］、均方誤差法［9］等。經(jīng)驗法需要依靠工程師的經(jīng)驗進行選取，具有很大的局限性。灰色關(guān)聯(lián)度法是根據(jù)選擇不同閾值所得超出量擬合曲線與樣本的分布曲線的接近程度的研究方法，具有較高的精確性；均方誤差法是運用統(tǒng)計學抽樣來模擬樣本中的極值，通過計算均方誤差進行選擇的方法。本文采用均方誤差法進行閾值選取。

本文中提出了一種基于廣義Pareto分布（generalized pareto distribution，GPD）函數(shù)的扭轉(zhuǎn)載荷時域外推方法。本文介紹廣義Pareto模型（Pareto模型為一種數(shù)學分布模型，即極值樣本服從Pareto分布）的建立過程，同時以商用車駕駛室穩(wěn)定桿的扭轉(zhuǎn)隨機載荷譜為外推對象，對比分析了外推前后載荷譜的穿級計數(shù)、功率譜密度、雨流圖和潛在偽損傷量。結(jié)果表明：所構(gòu)建的載荷外推算法能夠準確地實現(xiàn)扭轉(zhuǎn)載荷的外推。

1 基于GPD的扭轉(zhuǎn)載荷外推理論

扭轉(zhuǎn)載荷是部件發(fā)生疲勞失效的一種重要載荷形式。汽車結(jié)構(gòu)中的典型部件（如穩(wěn)定桿、傳動軸等）在服役過程中，在扭轉(zhuǎn)載荷作用下，其內(nèi)部將呈現(xiàn)出交變切向應(yīng)力，在較大幅值的應(yīng)力長時間作用后，結(jié)構(gòu)發(fā)生疲勞失效。與受水平載荷的典型部件不同，受純扭轉(zhuǎn)載荷的部件更多出現(xiàn)的是剪切疲勞失效，而非拉壓疲勞失效。因此，部件全壽命周期內(nèi)的扭轉(zhuǎn)隨機疲勞載荷譜是研究結(jié)構(gòu)疲勞性能的關(guān)鍵?？筛鶕?jù)采集的扭轉(zhuǎn)隨機載荷的統(tǒng)計規(guī)律，對未知載荷譜進行預(yù)測，實現(xiàn)目標里程下載荷譜的獲取，完成部件耐久性評估。

1.1 GPD函數(shù)及其性質(zhì)

對部件所受到的隨機載荷譜進行外推前，需要先進行載荷譜的預(yù)處理，如載荷的峰谷值提取。對于載荷譜中峰谷值超過某一閾值的部分，定義為載荷譜的極值（也叫超出量），見圖1。研究表明：隨機載荷譜的極值服從GPD分布［3］。廣義Pareto分布函數(shù)表達式為［10］：

式中：u為位置參數(shù)（實為所設(shè)置的閾值）；σ為尺度參數(shù)；ξ為形狀參數(shù)。其概率密度函數(shù)為

根據(jù)圖1，Xj（j＝1，2，3，…，N）表示提取峰谷值后的采樣點，N為載荷譜峰谷值的采樣點總數(shù)；設(shè)定閾值為u，其中采樣點中超過閾值的采樣點定義為極值樣本，記作Yk＝Xj－u（j＝1，2，3，…，N；k為樣本序列號，j為采集點序列號）。將載荷譜中的極值樣本分為上極值樣本和下極值樣本，即在載荷譜中提取出載荷譜的峰谷值，然后設(shè)定上閾值umax和下閾值umin，峰谷值中超過上閾值的部分就是載荷譜的上極值樣本，超過下閾值的部分就是載荷譜的下極值樣本，上下極值樣本可分別表示為

式中：為上極值樣本：為下極值樣本。

1.2 GPD函數(shù)閾值

基于廣義Pareto分布的載荷譜外推，閾值的選擇決定了廣義Pareto分布函數(shù)的參數(shù)，進而影響載荷譜外推結(jié)果。載荷譜閾值的確定首先需要基于極值樣本的均值超出函數(shù)確定一個取值區(qū)間，然后以形狀參數(shù)最小的均方誤差為目標，確定最佳閾值［9］。

1.2.1 經(jīng)驗均值的超出函數(shù)

通過經(jīng)驗均值超出函數(shù)初步確定一個閾值區(qū)間范圍，以形狀參數(shù)最小的均方誤差為目標，在該區(qū)間內(nèi)通過增量步迭代，實現(xiàn)最優(yōu)閾值的確定。經(jīng)驗均值超出函數(shù)為［4］

式中：Xi是極值樣本；Nu是樣本中的數(shù)量。根據(jù)GPD函數(shù)性質(zhì)，超出量均值函數(shù)為

閾值區(qū)間的確定，需要繪制經(jīng)驗均值超出函數(shù)en（u）的散點圖，隨后根據(jù)式（6）可知，當形狀參數(shù)和尺度參數(shù)一定時，經(jīng)驗均值超出函數(shù)和閾值具有線性關(guān)系。故在經(jīng)驗均值超出函數(shù)圖中找到呈線性關(guān)系的區(qū)間，則該區(qū)間為閾值范圍區(qū)間。

1.2.2 GPD函數(shù)最優(yōu)閾值的確定

基于GPD函數(shù)的極值樣本的矩分析，可以得到Pareto分布函數(shù)中形狀參數(shù)與矩之間的關(guān)系表達式，如下：

由式（7）可知，GPD函數(shù)的形狀參數(shù)ξ與所選定的閾值u相關(guān)。當閾值u偏大時，極值樣本的數(shù)量將會過少，導(dǎo)致樣本分布不穩(wěn)定；當閾值u偏小時，極值樣本的數(shù)量過多，將導(dǎo)致所得結(jié)果極大偏離估計值和真實值。因此，采用均方誤差MSE（Mean square error）描述閾值u的選擇情況。均方誤差綜合考慮閾值偏大或者偏小所帶來的誤差。均方誤差通過如下公式獲?。?］。

式中：Bia2（ξ）是形狀參數(shù)偏差；Var（ξ）是形狀參數(shù)方差。

形狀參數(shù)偏差用于描述估計值與實際值之前的巨大偏差，通過如下公式計算。

式中：b為抽樣次數(shù)。形狀參數(shù)的方差用于描述樣本分布的穩(wěn)定性，通過如下公式計算。

采用均方誤差估計閾值u時，使用自抽樣方法對極值樣本進行多次重復(fù)抽樣［11］，得到自抽樣法樣本組{B}i，進而得到多個相對應(yīng)的形狀參數(shù)估計值。其中，最小形狀參數(shù)均方誤差估計值所對應(yīng)的閾值為最佳閾值。

綜上，GPD函數(shù)最優(yōu)閾值的確定流程可總結(jié)如下（見圖2）：

圖2 閾值確定流程

1）在經(jīng)驗超出量函數(shù)確定的閾值區(qū)間中選擇一個閾值作為初始閾值u0，得到極值樣本，并記錄極值樣本中的數(shù)量Nj。

2）計算出形狀參數(shù)的初始估計值，記為ξ0。

3）使用自抽樣法對極值樣本隨機可重復(fù)選擇Nj個樣本數(shù)據(jù)，得到一組自抽樣樣本Bi。

4）計算出自抽樣樣本Bi對應(yīng)的形狀參數(shù)估計值ξi。

5）重復(fù)步驟3）和4），得到b組自抽樣樣本所對應(yīng)的形狀參數(shù)估計值組｛ξi｝。

6）計算自抽樣樣本組對應(yīng)的Bia（ξ）和Var（ξ）的估計量，然后運用式（8）計算MSE（ξi）。

7）重置初始閾值uj＝uj＋j·Δu（j＝1，2，…N），重復(fù)步驟1）—步驟6），計算得到閾值樣本組｛uj｝和對應(yīng)的均方誤差估計值組｛MSE（ξi）｝。

8）均方誤差估計值組｛MSE（ξi）｝中的最小值對應(yīng)的閾值uj為最佳閾值。

步驟5）中的抽樣次數(shù)b對閾值的分析結(jié)果有顯著影響。當抽樣次數(shù)太少，計算結(jié)果將不穩(wěn)定；而抽樣次數(shù)超過一定量時，計算結(jié)果將趨于穩(wěn)定。因此，要保證閾值的精確度，需要達到一定抽樣次數(shù)。

1.3 GPD函數(shù)形狀參數(shù)和尺度參數(shù)估計

GPD函數(shù)的形狀參數(shù)和尺度參數(shù)的常用方法有：矩估計法［12］、極大似然估計法［13］、概率加權(quán)矩法［14］等。其中，基于極大似然的估計法，首先對含有未知參數(shù)的分布函數(shù)建立對數(shù)似然函數(shù)，對似然函數(shù)求導(dǎo)并等于零，求解方程得到估計參數(shù)的估計量。具體而言，對極值樣本建立對數(shù)似然函數(shù)：

基于式（11）對形狀參數(shù)和尺度參數(shù)求偏導(dǎo)，得到關(guān)于形狀參數(shù)和尺度參數(shù)的方程組。

求解以上方程組，便可獲取形狀參數(shù)和尺度參數(shù)的估計值。

1.4 基于GPD的載荷譜外推流程

基于GPD函數(shù)的載荷譜外推的步驟總結(jié)如下：

1）輸入采集載荷譜并提取載荷譜峰谷值Xj。

2）選擇合適的上閾值umax和下閾值umin。

3）提取上閾值超出量Ymaxk和下閾值超出量Ymink。

4）采用GPD函數(shù)對超出量Yk進行擬合，獲得GPD函數(shù)的參數(shù)。

5）將擬合得到的參數(shù)還原到GPD函數(shù)中，然后通過該函數(shù)隨機生成數(shù)量相等的超出量。

6）生成的超出量替換原始載荷譜中的超出部分，從而獲得外推信號。

2 載荷譜采集與外推

2.1 載荷譜采集

商用車駕駛室在轉(zhuǎn)向等受到側(cè)向力的工況時，駕駛室的穩(wěn)定桿能夠防止車身發(fā)生過大橫向側(cè)傾，使駕駛室保持平衡。故駕駛室穩(wěn)定桿的本質(zhì)是一個橫置的扭桿彈簧，在功能上可以看成是一種特殊的彈性元件，其服役過程中主要受到扭轉(zhuǎn)載荷的作用。以商用車駕駛室穩(wěn)定桿為對象，首先介紹其扭轉(zhuǎn)載荷采集的方法，然后進行扭轉(zhuǎn)載荷外推研究與驗證。

2.2 穩(wěn)定桿載荷譜采集

在進行駕駛室穩(wěn)定桿扭轉(zhuǎn)載荷采集之前，需要分別進行貼片、組橋和標定工作，見圖3。

圖3 駕駛室穩(wěn)定桿扭轉(zhuǎn)載荷采集應(yīng)變片粘貼、組橋與標定

將貼片、標定好后的駕駛室穩(wěn)定桿裝車，在典型的駕駛工況上進行載荷譜采集，得到100 km的載荷譜數(shù)據(jù)，見圖4。

圖4 用戶典型工況100 km采集載荷譜數(shù)據(jù)

為了對外推算法進行驗證，將采集的100 km用戶城市工況載荷譜數(shù)據(jù)分為兩部分，一部分載荷譜的里程長度為50 km，用于基于廣義Pareto分布的扭轉(zhuǎn)載荷外推；另一部分載荷譜的里程長度為50 km，用于與外推扭轉(zhuǎn)載荷進行對比，驗證載荷外推的效果。根據(jù)圖4，駕駛室穩(wěn)定桿前50 km載荷數(shù)據(jù)中最大載荷為960.4 N·m，最小載荷為－1 378 N·m；駕駛室穩(wěn)定桿后50 km載荷數(shù)據(jù)中最大載荷為1 378 N·m，最小載荷為－1 192 N·m。

2.3 GPD函數(shù)參數(shù)擬合

2.3.1 閾值確定

根據(jù)式（8），采用均方誤差MSE來確定閾值u，但首先需要經(jīng)驗均值的超出函數(shù)初步確定一個閾值區(qū)間范圍，見圖5。

圖5 上下閾值的超出量分布圖

根據(jù)試驗載荷樣本，初步確定上閾值的區(qū)間范圍為350～400 N·m，下閾值的區(qū)間范圍為－330～－380 N·m，見圖6。初步確定上下閾值區(qū)間后，使用MSE來描述上下閾值區(qū)間內(nèi)的閾值選擇情況，上下閾值對應(yīng)的MSE分別見圖6（a）與圖6（b）。

圖6 上下閾值及其對應(yīng)均方誤差曲線

根據(jù)圖6，在指定的閾值區(qū)間內(nèi)，每個閾值所對應(yīng)的MSE隨閾值的增大而整體呈現(xiàn)波動上升的趨勢，但分別在351 N·m和－331 N·m時取得最小值。因此，最終確定上閾值為351 N·m，下閾值為－331 N·m。

2.3.2 參數(shù)擬合及評價

基于確定的上、下閾值，從駕駛室穩(wěn)定桿扭轉(zhuǎn)載荷數(shù)據(jù)樣本中，分別提取上、下閾值極值樣本并進行擬合，見圖7。

圖7 峰谷值極值樣本GPD函數(shù)擬合

根據(jù)圖7，GPD函數(shù)與峰谷值極值樣本擬合時，較小出現(xiàn)偏離參考線的情況，說明GPD函數(shù)擬合情況較為精確，具有較小的擬合參數(shù)偏差。為了對擬合結(jié)果進行量化評價，采用擬合優(yōu)度對擬合結(jié)果進行評判［15］。

擬合優(yōu)度是指回歸曲線對離散數(shù)據(jù)的擬合程度，判定系數(shù)R2度量回歸曲線與離散數(shù)據(jù)之間的吻合程度。R2的取值范圍為［0，1］，當R2值越接近1，說明回歸曲線對離散數(shù)據(jù)的擬合程度越好；反之，說明回歸曲線對離散數(shù)據(jù)的擬合程度越差。R2值為0.7～1，相關(guān)性強；R2值為0.5～0.7，相關(guān)性一般；R2值為0～0.5，相關(guān)性弱［16］。

表1 基于峰谷值超出量的GPD函數(shù)擬合參數(shù)

根據(jù)表1，峰值極值樣本的相關(guān)系數(shù)R2值為0.875 2，谷值極值樣本的相關(guān)系數(shù)R2值為0.793，說明通過擬合得到峰谷值超出量的形狀參數(shù)和尺度參數(shù)，與極值樣本具有強相關(guān)性，能夠保留極值樣本的統(tǒng)計學特性。

3 載荷譜外推結(jié)果與分析

3.1 載荷譜外推

根據(jù)載荷譜外推步驟，通過GPD函數(shù)模型隨機生成數(shù)量相等的峰值和谷值超出量，并替換原始載荷譜中的相應(yīng)峰谷值，獲得外推載荷譜。對駕駛室穩(wěn)定桿前50 km扭轉(zhuǎn)載荷進行外推，得到的外推扭轉(zhuǎn)載荷見圖8。

圖8 基于GPD函數(shù)的外推載荷譜

根據(jù)圖8，基于GPD函數(shù)外推得到的載荷譜，其極值的絕對值與原始載荷譜相比變大。極大值從原始數(shù)值960.4 N·m增大為1 939.6 N·m，極小值從原始數(shù)值－1 378 N·m增大為－1 409.4 N·m。

3.2 載荷譜外推結(jié)果分析

為進一步對比扭轉(zhuǎn)載荷的外推效果，從穿級計數(shù)、雨流計數(shù)、功率譜密度以及偽損傷分析等方面，對外推前后載荷譜進行對比，驗證所提出的外推算法。

3.2.1 穿級計數(shù)分析

對前50 km采集到的原始載荷譜、GPD函數(shù)外推得到的載荷譜，以及后50 km采集到的原始載荷譜進行穿級計數(shù)對比分析，觀察其幅值和頻次的變化情況，見圖9。

圖9 原始載荷譜與外推載荷譜穿級計數(shù)曲線

根據(jù)圖9，前50 km采集的原始載荷譜在基于GPD函數(shù)外推后，所得到的外推載荷譜在幅值較小處的頻次基本吻合，在幅值較大處的頻次有所增加。主要原因：原始載荷譜外推后其極值的絕對值大于原始載荷譜極值的絕對值，由于外推載荷極值大小生成具有隨機性，是生成少量較大極值導(dǎo)致。

基于GPD函數(shù)外推得到的載荷譜，在負值扭轉(zhuǎn)載荷中與后50 km采集的原始載荷譜更加吻合。但在正值扭轉(zhuǎn)載荷中，外推載荷譜與實際采集的載荷譜之間具有一定的偏差，主要體現(xiàn)在：外推載荷譜出現(xiàn)了少量較大幅值的循環(huán)，而采集的后50 km原始載荷譜中間幅值的循環(huán)較多。

3.2.2 雨流分析

對駕駛室穩(wěn)定桿前50 km 原始載荷譜、前50 km外推載荷譜和后50 km原始載荷譜進行雨流計數(shù)，繪制出From-to雨流圖，見圖10。

圖10 原始載荷譜與外推載荷譜雨流計數(shù)

根據(jù)圖10（a）和（b），前50 km采集的原始載荷譜相比后50 km采集的原始載荷譜，大幅值信號明顯較少。圖10（c）的外推載荷譜雨流圖，有一個大幅值循環(huán)出現(xiàn)，表明基于GPD函數(shù)的外推得到的外推載荷譜有較大幅值的出現(xiàn)?？傮w來看，前50 km GPD函數(shù)外推得到的外推載荷譜與后50 km采集到的原始載荷譜的效果更加一致。

3.2.3 功率譜密度分析

將前50 km外推載荷譜與后50 km原始載荷譜進行能量變化分析，將其功率譜密度進行對比。

根據(jù)圖11，在1 Hz處出現(xiàn)峰值，由振動理論可知，該頻率為駕駛室穩(wěn)定桿的固有頻率。在20 Hz之前，前50 km外推載荷譜能量與后50 km原始載荷譜能量基本吻合；但20 Hz之后，前50 km外推載荷譜能量與后50 km原始載荷譜能量有一定差異，主要原因在于外推后的載荷譜中增加了部分大幅值信號，而低頻對應(yīng)小幅值，高頻對應(yīng)大幅值，因此使得外推后的載荷譜功率譜密度的高頻部分高于原始載荷譜功率譜密度。但外推前后載荷譜的能量主要集中在前20 Hz，因此外推后的載荷譜能量與原始載荷譜基本一致。

圖11 外推前后功率譜密度變化曲線

3.2.4 偽損傷分析

為了進一步驗證載荷外推方法的準確性，對前50 km原始載荷譜、前50 km外推載荷譜和后50 km原始載荷譜進行偽損傷分析。運用ncode軟件進行偽損傷分析，商用車扭桿的材料相關(guān)參數(shù)見表2。

表2 商用車扭桿的材料屬性

分析結(jié)果如表3所示，前50 km外推載荷譜潛在偽損傷量略大于前50 km原始載荷譜，原因是外推后載荷譜的載荷循環(huán)次數(shù)略大于原始載荷譜；后50 km原始載荷譜潛在偽損傷量小于前50 km外推載荷譜，但循環(huán)次數(shù)大于前50 km原始載荷譜，其原因為后50 km原始載荷循環(huán)在小應(yīng)力幅處的循環(huán)次數(shù)高于前50 km外推載荷，在大應(yīng)力幅值處的循環(huán)次數(shù)低于前50 km外推載荷。此外，前50 km、后50 km原始載荷譜的循環(huán)次數(shù)與路面狀況的隨機性存在一定關(guān)系。后50 km原始載荷譜與外推載荷譜潛在偽損傷量差異性較小，表明GPD函數(shù)載荷外推在誤差允許范圍內(nèi)，具有較好的結(jié)果。

表3 3種載荷潛在偽損傷量計算結(jié)果

4 結(jié)論

通過所構(gòu)建的載荷外推算法進行外推，其外推載荷譜極值的絕對值均大于原始載荷譜極值的絕對值。外推后的載荷譜中增加了部分大幅值信號，而低頻對應(yīng)小幅值，高頻對應(yīng)大幅值，因此使得外推后的載荷譜功率譜密度的高頻部分高于原始載荷譜功率譜密度，外推載荷潛在偽損傷量與原始載荷潛在偽損傷量差異較小。故基于GPD函數(shù)外推算法，所得外推結(jié)果與原始載荷譜基本一致，GPD函數(shù)能夠準確實現(xiàn)商用車駕駛室穩(wěn)定桿扭轉(zhuǎn)載荷外推。因此，本文外推算法對商用車作業(yè)時駕駛室穩(wěn)定桿承受扭轉(zhuǎn)載荷具有較好適應(yīng)性。