基于核密度估計的拖拉機傳動軸載荷外推方法

李淑艷 翟友邦 王小龍 李若晨 楊子涵 宋正河

(中國農業大學 工學院，北京 100083)

載荷譜是機械零部件壽命設計和可靠性分析的基礎[1]。傳動軸作為拖拉機傳動系的主要組成部分，承擔著從變速箱向驅動橋傳遞動力的功能，是影響拖拉機整機可靠性的關鍵零部件。編制傳動軸載荷譜，可以提高疲勞試驗的可信度，對提高傳動系可靠性具有重要意義。出于成本考慮，載荷測試試驗時間相對較短，需要對測試載荷進行外推獲得全壽命載荷譜。

常用的外推方法有時域外推和雨流矩陣外推[2]。其中，時域外推能夠保留載荷時間順序，適用于平穩性載荷，當載荷表現出非平穩性時，其計算量較大。雨流矩陣外推雖不能保留載荷時間順序，但其計算量小，對平穩和非平穩載荷均適用，具有較強工程應用價值。由于拖拉機作業環境復雜，導致傳動軸承受非平穩性載荷，因此適合采用雨流矩陣外推的方法。

雨流矩陣外推常用的方法有參數法和非參數法。參數法通常采用正態分布擬合均值頻次，威布爾分布擬合幅值頻次。賈海波[3]基于單參數法對輪式裝載機傳動系載荷進行外推，在此基礎上獲得了均幅值二維載荷譜，并給出了疲勞試驗一維程序譜加載方案。已有研究[4-5]通過參數法對輪式裝載機半軸載荷進行外推，其中，翟新婷提出的混合分布的載荷譜編制方法對載荷多峰分布的擬合難題拓寬了解決思路，但擬合精度仍需提高。陳東升等[6]以軍用車輛傳動系統為研究對象，采用雙參數雨流計數法，對載荷的幅值和均值進行計數，并將測試載荷外推得到全壽命歷程二維設計載荷譜，建立了可用于疲勞壽命估算的八級程序載荷譜。采用基于參數法的雨流矩陣外推方法時，由于非平穩性載荷具有單峰和多峰的特點，存在擬合不通過或分布函數難以確定的局限。非參數法是基于核密度估計(Kernel density estimation，KDE)的思想對雨流矩陣進行外推，可以直接得出任意載荷的密度值，克服分布函數擬合的缺陷，為多峰復雜載荷的精確擬合提供了解決辦法。李鶯鶯等[7]采用非參數估計方法對挖掘機液壓泵載荷進行外推，得到全生命周期內每個載荷循環可能出現的頻次，而且可以保證每個遲滯回環的結構不被破壞。張曉晨等[8]選擇基于載荷擴展的非參數外推方法分別對液壓挖掘機的4個作業段進行外推，獲得全壽命下的長期載荷譜。李研[9]提出了基于全帶寬非參數外推模型和載荷擴展的非參數雨流外推方法，針對輪式裝載機，根據實測載荷數據進行了試驗驗證與載荷譜編制工作，并指出帶寬對非參數外推精度的局限性。Johannesson.P等[10]將極值理論與非參數法結合起來，對汽車半軸載荷進行了外推處理。

雖然國內外都對非參數外推開展了相關研究，但外推精度依然亟待提高。針對上述問題，本研究擬以拖拉機傳動軸為研究對象，搭建傳動軸無線扭矩測試系統并進行載荷測試，分析測試載荷分布特征，運用拇指法(Rule-of-thumb，ROT)求帶寬，提出一種改進的基于核密度估計的載荷外推方法對測試載荷進行外推，并從參數統計、擬合度和偽損傷3個方面對所提出的方法的精度和有效性進行驗證。

1 拖拉機傳動軸田間試驗

試驗樣機為LX2204拖拉機，所搭載的犁耕機型號為MULTI-MASTER153T。在傳動軸中心附近粘貼應變片，同時無線扭矩節點固定在傳動軸上以確保傳感器與傳動軸的相對位置保持不變。選用北京必創科技有限公司生產的TQ201型無線扭矩采集器，其使用應變片為BF350-3BA型，基于TQ201搭建的拖拉機傳動軸無線扭矩測試系統見圖1。

圖1 拖拉機傳動軸無線扭矩測試系統Fig.1 Tractor propeller shaft torque test system

田間試驗地點為河南省洛陽市孟津縣金村面積約13 hm2的砂土農田，玉米秸稈留茬高度約10 cm。采集了犁耕、旋耕、運輸和聯合耕整等工況的傳動軸載荷試驗數據，采樣頻率為100 Hz。為保證測試載荷能真實反映作業工況，根據GB/T 14225—2008《鏵式犁》中的規定，試驗過程由駕駛員按其習慣進行操作。試驗中拖拉機的檔位為中三檔，行駛速度為5～10 km/h；耕深和作業寬幅分別為0.32和2.1 m；碎土率為99.2%。對田間作業環境進行測量，風速為4.6 m/s，環境濕度和溫度分別為21%和25.6 ℃；距土壤表面5、10和15 cm處的含水率分別為6.28%、9.84%和12.2%。

2 基于核密度估計的載荷外推方法

本研究采用雨流矩陣外推中基于核密度估計理論的非參數法和按比例外推方法對測試載荷進行外推處理，并對非參數法在外推過程中帶寬的選取進行探究，載荷的外推流程見圖2。

圖2 基于核密度估計的載荷外推方法流程Fig.2 Process of load extrapolation method based on kernel density estimation

2.1 信號預處理

田間試驗共獲得8組犁耕工況傳動軸載荷數據，受測試環境和測試系統等原因影響，田間測試信號中不可避免地會存在干擾噪聲。通過多項式擬合去除測試信號中的趨勢項、幅值門限法和動態標準差法去除信號中的奇異點以及巴特沃斯低通濾波器濾除高頻載荷，由于高頻載荷大部分是小載荷，因此低通濾波可以將對結構損傷貢獻很小的小載荷濾除。圖3為預處理后的部分犁耕工況傳動軸測試載荷時間歷程。

圖3 犁耕工況傳動軸扭矩時間歷程Fig.3 Torque time history of propeller shaft in ploughing conditions

2.2 核密度估計

采用核密度估計理論對雨流矩陣進行外推屬于非參數法。核密度估計能較準確估計任意分布載荷數據的密度值，當帶寬選擇恰當時，密度值能達到較高的精度，核密度估計的具體步驟為[11]：

1)采用雨流計數法將載荷-時間歷程轉換成雨流矩陣。

2)確定核函數K(x,y)和帶寬h，構建核密度估計模型f(x,y)。

3)根據核密度估計模型和敏感系數計算自適應因子λi：

(1)

式中：α為敏感系數，經驗值為0.5；n為樣本量，xi、yi分別為雨流矩陣初值(From)和終值(To)。

4)根據自適應因子確定自適應帶寬的核密度估計模型：

(2)

首先采用核密度估計模型對雨流矩陣進行密度估計獲得密度矩陣，然后采用Monte Carlo法對密度矩陣進行隨機數模擬，生成外推載荷。

2.3 帶寬確定

從式(2)可以看出，核函數和帶寬影響核密度估計精度。研究表明[12-13]，核函數對核密度估計精度影響較小，帶寬是影響核密度估計的主要因素。帶寬過大，會導致載荷特性丟失；帶寬過小，則會在密度中引入大量虛假波型，造成密度估計不穩定。本研究結合拖拉機傳動軸載荷特征確定ROT帶寬計算方法。

2.3.1R0T帶寬計算方法

本研究選取高斯核函數作為核密度估計模型的核函數，并在核密度估計模型的基礎上確定帶寬計算方法。針對傳動軸測試載荷，采用四點雨流計數法將測試載荷轉換成初值(From)- 終值(To)雨流矩陣。雨流矩陣中存在單峰載荷和多峰載荷(圖4)。ROT能根據單峰載荷和多峰載荷的特性，通過載荷樣本的標準差和標準四分位距(Standardized interquartile range)確定出最優帶寬[14-16]。

圖4 測試載荷雨流矩陣分布特征Fig.4 Test load rain flow matrix distribution characteristics

基于高斯核函數的核密度估計模型的ROT求帶寬h的公式為：

(3)

式中：A為高斯分布的最小值；σ為樣本標準差；IQR為數據的標準四分位距。

拖拉機傳動軸載荷具有大載荷較為分散，中小載荷相對集中的特點。對犁耕工況傳動軸測試載荷雨流矩陣中的載荷進行簡化(圖5)。當用式(3)計算的帶寬應用于a點的密度估計時，帶寬會偏小。在以帶寬為半徑的圓內計算a點的概率密度分布，并進行Monte Carlo隨機數模擬。如圖6所示，由于帶寬偏小，產生的隨機數會過于集中?？梢赃m當增大帶寬以提高隨機數的分散性，通過減小樣本量來增大帶寬。參照時域外推選取大載荷閾值的原則和實際外推效果，選擇總樣本量的80%計算帶寬，則修正樣本量N為：

a代表分散性較強的大載荷，b代表分布相對集中的中小載荷。Point a represents a large load with strong dispersion, point b represents a relatively concentrated medium and small load.圖5 犁耕工況測試載荷雨流矩陣分布點Fig.5 Rain flow matrix distribution points of test load in plowing conditions

a*代表對a點的大載荷進行密度估計。a* represents the density estimation of the large load at point a.圖6 大載荷核密度估計Fig.6 Rain flow matrix large load kernel density estimation

N=[0.8×n]+1

(4)

2.3.2核密度估計應用范圍

應用式(1)計算出的載荷分布集中的區域，其自適應因子較小，即發生概率較大，針對這類載荷采用按比例外推法可以保留其原分布，且外推后仍在該區域，符合原載荷較集中的特征。同時采用線性外推后，在Monte Carlo模擬時略去此區域，計算速度可得到大幅提升。

對于雨流矩陣分散性較強(大載荷)的區域使用基于KDE的非參數法外推。在Monte Carlo模擬時會在其周圍產生新的分布點(圖7)，新分布點符合原載荷分散性強的特性。

a’ 代表Monte Carlo模擬產生的新分布點。a’ represents the new distribution point generated by Monte Carlo simulation.圖7 KDE外推結果Fig.7 The results of KDE extrapolation

樣本量也是影響核密度估計精度的一個因素，為確保有較高的外推精度，需保證非參數法有足夠多的樣本量，因此需要在載荷的分散性和樣本量之間尋求一個合適的閾值。根據拖拉機傳動軸在各工況下的載荷分布特征及載荷對疲勞損傷的貢獻度確定閾值。

1)計算雨流矩陣中載荷的幅值a：

(5)

式中：F為雨流矩陣的初值；T為雨流矩陣的終值。

2)定義非參數法和按比例外推法之間的閾值t：

t=[0.1max(a)]+1

(6)

通過上述分析和修正，用于拖拉機傳動軸載荷外推的方法為：

(7)

式中：k為外推倍數；Ns為小載荷循環數；G(x,y)為采用基于KDE的非參數法的載荷外推函數。

3 載荷外推

以圖3所示犁耕工況為例進行載荷外推計算。因測試載荷為扭矩，雨流計數是基于材料應力遲滯回線原理形成的，故將扭矩轉換成應力進行計數處理，轉換公式為：

(8)

式中：τ為應力，Pa；T為扭矩，N·m；Wt為抗扭截面系數，Wt=πd3/16，d為軸的直徑，m。

將犁耕工況傳動軸扭矩轉換成應力(圖8)，提取測試應力數據的轉折點并用四點雨流計數法得到雨流載荷循環的初值、終值及其頻次信息，對載荷進行分級并編制成雨流矩陣。本研究將雨流載荷分成64級，編制成64×64雨流矩陣見圖9。雨流矩陣的原始樣本量n=1 488，根據式(4)計算出修正樣本量N=1 191。根據式(5)和式(6)進一步計算出閾值t=8.2 MPa，結合式(3)可計算出最優帶寬為2.709和2.691。非參數法的外推參數見表1。

表1 基于核密度估計的非參數法外推參數Table 1 The extrapolation of parameters based on nonparametric method of KDE

圖8 梨耕工況傳動軸應力時間歷程Fig.8 Stress time history of propeller shaft in ploughing conditions

圖9 基于應力數據編制的雨流矩陣Fig.9 Rainflow matrix based on stress data

根據閾值t劃分非參數法外推載荷和按比例外推載荷。非參數法外推載荷見圖10(a)，樣本數量為318，按比例外推載荷見圖10(b)，樣本數量為1 170。

圖10 犁耕工況傳動軸載荷外推結果Fig.10 The result of different extrapolation methods of propeller shaft load in ploughing condition

取敏感系數α=0.2，采用上述核密度估計模型對雨流矩陣進行密度估計，并用Monte Carlo模擬產生與樣本量相同的載荷。把非參數法外推的載荷和按比例外推的載荷合并形成一個雨流矩陣見圖11,對比圖9可以看出，外推1倍的雨流矩陣載荷分布與原始雨流矩陣載荷分布具有較高的相似度，保留了原有載荷的分布特征。相比于整段載荷使用核密度估計的傳統非參數法，外推的運算速度得到大幅提高，當外推至全壽命頻次(106)時，運算時間從0.502 310 s縮短至0.230 247 s，速度提高了54.16%。

圖11 犁耕工況傳動軸總體外推結果Fig.11 Total result of extrapolation of propeller shaft in ploughing condition

4 外推結果驗證

當外推方法選擇合理時，外推載荷與測試載荷應具有相同的特征，本研究分別從統計參數、擬合度和偽損傷 3 個方面檢驗外推載荷的有效性。

4.1 統計參數檢驗

選取統計參數均值、標準差和最大值作為傳動軸測試載荷與外推載荷的評判依據。對犁耕工況，采用本研究提出的基于核密度估計的載荷外推方法得到的外推載荷，其均值、標準差和最大值與測試載荷的誤差分別為2.2%、0.87%和6.3%，均在工程應用允許誤差范圍內。旋耕工況、運輸工況和聯合耕整工況的均值、標準差和最大值誤差這3種統計參數均小于10%。綜上表明外推載荷能較好地保留測試載荷的統計特性。各工況測試載荷與外推載荷的統計參數見表2。

表2 不同工況下傳動軸載荷統計參數Table 2 Statistical parameters of propeller shaft load under different working conditions MPa

4.2 擬合度檢驗

采用擬合度校驗外推結果。載荷的循環次數和幅值是疲勞損傷的主要來源，通過繪制幅值累計頻次曲線可以反映載荷的幅值和循環次數的統計結果。由外推載荷和測試載荷繪制成的雨流載荷幅值累計頻次對比曲線見圖12。

圖12 犁耕工況幅值累計頻次曲線Fig.12 Accumulated frequency curve of amplitude of ploughing condition

為判斷外推載荷與測試載荷之間的差異，對曲線擬合度進行量化，用決定系數R2(Coefficient of determination)計算曲線之間的擬合度。對犁耕工況，外推載荷幅值的R2達到0.998，其中前100個大幅值載荷的R2為0.978，表明該方法保留了拖拉機傳動軸載荷的分布特性。采用旋耕工況、運輸工況和聯合耕整工況下傳動軸測試載荷進一步驗證本研究所提外推方法的合理性，并與傳統非參數法進行對比，結果表明該外推方法對傳動軸不同工況載荷具有較高的適應性，外推的整段載荷和前100個大幅值載荷的R2均大于0.9，與測試載荷高度擬合。與傳統非參數法相比，由本研究所提外推方法得到的R2均有較大幅度的提高，外推效果明顯優于傳統非參數法。決定系數R2見表3。

表3 不同工況下外推載荷的決定系數Table 3 Determination coefficient R2 of extrapolated load under different working conditions

4.3 偽損傷檢驗

外推載荷的最終目的是編制疲勞試驗加載譜。在機械工程領域，分析復雜的交變載荷時，常通過計數法和損傷累積來計算疲勞。采用S-N曲線疲勞壽命模型計算構件的疲勞壽命，循環壽命與恒幅值應力的關系[17]為：

(9)

式中：α為材料屬性參數；β為損傷指數；Si為恒幅值應力；Ni為循環壽命。

工程上常用Palmgren-Miner線性損傷理論計算偽損傷，該理論是基于S-N曲線各應力循環次數對材料造成偽損傷的線性累積結果[18]，表達式為：

(10)

式中：Di為總偽損傷；ni為材料在Si載荷作用下的循環次數；Ni為載荷Si作用下的疲勞壽命循環次數。

為便于對偽損傷進行比較，定義偽損傷比為外推載荷對材料造成的偽損傷與測試載荷對材料造成的偽損傷的比值。偽損傷比越接近1，外推載荷與測試載荷的相似度越高?；赑almgren-Miner線性損傷疲勞累積理論在Ncode軟件中搭建偽損傷計算模型。

由本研究所提外推方法得到的犁耕工況外推載荷，其偽損傷比為0.999，與測試載荷具有高度一致性。同時計算旋耕工況、運輸工況和聯合耕整工況傳動軸載荷的偽損傷比，其數值分別為0.992、0.998 和0.996，均能較完整地保留測試載荷的損傷量。

5 結 論

本研究根據拖拉機傳動軸載荷特征，在采用ROT求得帶寬的基礎上，提出了一種基于核密度估計的載荷外推方法，并從統計參數、外推擬合度和偽損傷3個方面驗證了方法的有效性。主要結論如下：

1)與傳統非參數法相比，當載荷外推至全壽命頻次(106)時，采用本研究所提外推方法可將計算時間縮短54.16%，當測試載荷量較多時，能夠有效縮短外推時間，提高計算效率。

2)將犁耕工況的外推載荷與測試載荷進行對比分析。統計參數檢驗表明，均值、標準差和最大值的誤差分別為2.2%、0.87%和6.3%，外推載荷很好地保留了測試載荷的統計特征；擬合度檢驗表明，前100個大幅值載荷的R2為0.978，整段載荷R2達0.998，外推載荷很好地保留了測試載荷的分布特征；偽損傷檢驗表明，偽損傷比達到0.99，外推載荷很好地保留了測試載荷的損傷量。

3)用本研究所提外推方法對旋耕工況、運輸工況和聯合耕整工況傳動軸載荷進行外推，統計參數誤差均小于10%，R2均大于0.9，偽損傷比均大于0.99，所提方法適用于拖拉機傳動軸不同工況載荷外推。