極值服從廣義Pareto分布的扭轉載荷外推方法研究

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

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

0 引言

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

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

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

1 基于GPD的扭轉載荷外推理論

扭轉載荷是部件發生疲勞失效的一種重要載荷形式。汽車結構中的典型部件（如穩定桿、傳動軸等）在服役過程中，在扭轉載荷作用下，其內部將呈現出交變切向應力，在較大幅值的應力長時間作用后，結構發生疲勞失效。與受水平載荷的典型部件不同，受純扭轉載荷的部件更多出現的是剪切疲勞失效，而非拉壓疲勞失效。因此，部件全壽命周期內的扭轉隨機疲勞載荷譜是研究結構疲勞性能的關鍵?？筛鶕杉呐まD隨機載荷的統計規律，對未知載荷譜進行預測，實現目標里程下載荷譜的獲取，完成部件耐久性評估。

1.1 GPD函數及其性質

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

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

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

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

1.2 GPD函數閾值

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

1.2.1 經驗均值的超出函數

通過經驗均值超出函數初步確定一個閾值區間范圍，以形狀參數最小的均方誤差為目標，在該區間內通過增量步迭代，實現最優閾值的確定。經驗均值超出函數為［4］

式中：Xi是極值樣本；Nu是樣本中的數量。根據GPD函數性質，超出量均值函數為

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

1.2.2 GPD函數最優閾值的確定

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

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

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

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

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

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

綜上，GPD函數最優閾值的確定流程可總結如下（見圖2）：

圖2 閾值確定流程

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

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

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

4）計算出自抽樣樣本Bi對應的形狀參數估計值ξi。

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

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

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

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

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

1.3 GPD函數形狀參數和尺度參數估計

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

基于式（11）對形狀參數和尺度參數求偏導，得到關于形狀參數和尺度參數的方程組。

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

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

基于GPD函數的載荷譜外推的步驟總結如下：

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

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

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

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

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

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

2 載荷譜采集與外推

2.1 載荷譜采集

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

2.2 穩定桿載荷譜采集

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

圖3 駕駛室穩定桿扭轉載荷采集應變片粘貼、組橋與標定

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

圖4 用戶典型工況100 km采集載荷譜數據

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

2.3 GPD函數參數擬合

2.3.1 閾值確定

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

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

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

圖6 上下閾值及其對應均方誤差曲線

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

2.3.2 參數擬合及評價

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

圖7 峰谷值極值樣本GPD函數擬合

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

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

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

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

3 載荷譜外推結果與分析

3.1 載荷譜外推

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

圖8 基于GPD函數的外推載荷譜

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

3.2 載荷譜外推結果分析

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

3.2.1 穿級計數分析

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

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

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

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

3.2.2 雨流分析

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

圖10 原始載荷譜與外推載荷譜雨流計數

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

3.2.3 功率譜密度分析

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

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

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

3.2.4 偽損傷分析

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

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

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

表3 3種載荷潛在偽損傷量計算結果

4 結論

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