基于最優(yōu)分布擬合的拖拉機三點懸掛牽引力載荷譜編制與臺架試驗

王 玲，宗建華，王 禹，扶蘭蘭，毛 旭，王書茂

（中國農(nóng)業(yè)大學工學院，北京 100083）

0 引 言

拖拉機三點懸掛裝置是拖拉機關鍵零部件之一，各種不同類型農(nóng)機具通過三點懸掛裝置快速連接到拖拉機上，有效地提高了農(nóng)田耕作效率。由于田間作業(yè)工況的復雜多變，拖拉機懸掛點承受較大的隨機載荷，容易發(fā)生疲勞失效，直接影響拖拉機安全及作業(yè)效率。因此拖拉機三點懸掛裝置作業(yè)載荷特性與分布研究對于預測零件壽命、產(chǎn)品可靠性分析等具有重要意義。

基于實測數(shù)據(jù)編制的載荷譜是疲勞壽命分析和可靠性試驗的重要依據(jù)。目前，采用參數(shù)外推來預測不同工況下的全壽命載荷譜應用較為廣泛?；趨?shù)外推方法，賈海波等分別開展了輪式裝載機傳動系統(tǒng)、推土機車架縱梁、軍用車輛傳動系等零部件載荷譜測試及編制方法研究，翟友邦等以拖拉機動力輸出軸及懸掛裝置為研究對象，基于實測載荷數(shù)據(jù)進行了載荷譜與加速譜的編制。參數(shù)外推方法對單峰且分布均勻的載荷分布具有較好的擬合效果，對于多峰及概率分布不明確的載荷數(shù)據(jù)，傳統(tǒng)參數(shù)外推編制方法中單分布難以實現(xiàn)載荷多峰分布的擬合問題，混合分布的研究更接近實際作業(yè)工況的載荷分布，基于此，翟新婷等基于混合分布函數(shù)采用雙分布對測試載荷分布進行更精確的擬合，與單分布對比，雙分布擬合效果更好；田維波等將EM（Expectation-Maximization）算法分別用于混合高斯分布和混合威布爾分布的數(shù)據(jù)擬合，編制裝載機工作裝置及數(shù)控刀架載荷譜。以上研究多基于載荷分布確定的工況開展研究，對于實際作業(yè)工況中載荷隨機多變，單峰、多峰等分布難以確定的特征，上述研究無法實現(xiàn)各工況載荷分布的最優(yōu)擬合，直接影響外推極值精度與載荷編制效果。

本文以拖拉機三點懸掛牽引力為研究對象，在傳統(tǒng)參數(shù)外推基礎上提出了一種基于最優(yōu)分布擬合的載荷譜編制方法，通過建立混合分布模型并確定最優(yōu)基函數(shù)個數(shù)，采用EM算法進行最優(yōu)參數(shù)估計實現(xiàn)了極值頻次外推與載荷譜編制。基于此，開發(fā)拖拉機三點懸掛牽引力室內(nèi)臺架加載系統(tǒng)，實現(xiàn)外推載荷譜的高精度模擬加載與復現(xiàn)，擬為拖拉機關鍵零部件載荷測試及疲勞預測試驗提供理論依據(jù)和平臺支撐。

1 拖拉機田間試驗載荷采集

拖拉機三點懸掛牽引力田間試驗在農(nóng)業(yè)農(nóng)村部南京自動化機械研究所白馬基地開展，試驗樣機選用耕王RS-1254拖拉機配套1s-200深松機，在麥茬30 cm左右的麥田進行深松作業(yè)，試驗作業(yè)平均速度約為4.2 km/h。利用拖拉機無線采集系統(tǒng)對三點懸掛牽引力載荷進行采集，田間試驗及設備布置如圖1所示。

圖1 田間試驗及設備布置 Fig.1 Field experiment and equipment layout

拖拉機無線采集系統(tǒng)主要由軸銷力傳感器、數(shù)據(jù)采集箱、無線遠程傳輸網(wǎng)橋及上位機監(jiān)控軟件等組成。由于拖拉機作業(yè)環(huán)境惡劣、田間基站網(wǎng)絡信號差等問題，本文采用無線網(wǎng)橋?qū)崿F(xiàn)數(shù)據(jù)采集箱與上位機監(jiān)控軟件的遠程數(shù)據(jù)傳輸，穩(wěn)定傳輸距離可達10 km；三點懸掛牽引力載荷頻率為3～5 Hz，根據(jù)實際工程應用需求，采集頻率應是載荷頻率的3～4倍，因此數(shù)據(jù)采集箱采樣頻率為20 Hz。數(shù)據(jù)采集箱將數(shù)據(jù)無線傳輸至上位機監(jiān)控軟件，實現(xiàn)數(shù)據(jù)無線采集、顯示與保存等。

由于作業(yè)環(huán)境復雜，采集的三點懸掛牽引力信號可能會存在異常數(shù)據(jù)，直接影響信號分析以及載荷譜的外推精度。為了獲取三點懸掛牽引力載荷變化規(guī)律的真實信號，采用幅值門限、梯度門限、方差檢測相結(jié)合的方法對載荷數(shù)據(jù)去除奇異點并進行濾波，并對處理前后數(shù)據(jù)進行統(tǒng)計特性分析，預處理后數(shù)據(jù)載荷歷程如圖2所示，預處理前后載荷極值沒有變化，最大值和最小值分別為4.159 kN、0.259 kN，均值、標準差、方差分別降低0.4%、1.0%、1.7%。

圖2 預處理后三點懸掛牽引力載荷歷程 Fig.2 Traction load course of three-point suspension after pretreatment

2 基于最優(yōu)分布擬合的參數(shù)估計與驗證

混合分布模型是指總體樣本服從個來自不同分布特性的單樣本分布函數(shù)，其分布模型表達式為

混合分布模型具有普適性，但在實際作業(yè)工況中載荷分布無法確定，因此很難使用混合分布模型進行分布擬合。在參數(shù)外推均值服從正態(tài)分布、幅值服從威布爾分布基礎上，可以認為單峰分布或多峰分布的復雜載荷，其均值服從混合高斯分布，幅值服從混合威布爾分布，這樣既可以降低計算復雜程度，也可很好地涵蓋各工況分布特性載荷的擬合外推。

2.1 載荷統(tǒng)計計數(shù)與均幅值獨立性檢驗

對于隨機性較強的載荷，在進行載荷譜外推編制時，通常應用雨流計數(shù)法對測試載荷進行循環(huán)統(tǒng)計計數(shù)轉(zhuǎn)換成雨流域載荷，本文采用四點雨流計數(shù)法對處理后載荷數(shù)據(jù)進行循環(huán)載荷頻次統(tǒng)計，得到均值-幅值雨流矩陣以及均幅值頻次圖如圖3所示。

圖3 懸掛牽引力載荷雨流矩陣及均值幅值頻次 Fig.3 Suspension traction rainflow matrix and mean amplitude frequency

為獲得外推極值和頻次，需要根據(jù)雨流矩陣以及均值、幅值頻次數(shù)據(jù)，求出均值、幅值聯(lián)合概率分布函數(shù)。若均值、幅值相互獨立，其單獨概率密度函數(shù)相乘即可得到聯(lián)合概率密度。若均值、幅值不相互獨立，則需要通過其他方法來確定其聯(lián)合概率分布函數(shù)。

2.2 最優(yōu)基函數(shù)個數(shù)確定

混合分布基函數(shù)個數(shù)即載荷分布是混合分布參數(shù)估計的前提和關鍵，且基函數(shù)個數(shù)直接影響混合分布的擬合效果，基函數(shù)個數(shù)過少會造成欠擬合，無法完整體現(xiàn)載荷數(shù)據(jù)分布，個數(shù)過多形成過擬合，部分數(shù)據(jù)不符合載荷分布。目前常用的基函數(shù)個數(shù)確定方法是將基函數(shù)從1開始遞增，逐一進行混合分布參數(shù)估計，通過與原始數(shù)據(jù)擬合優(yōu)度對比，得到最優(yōu)基函數(shù)個數(shù)。通過總結(jié)，對于較為復雜的載荷分布，混合高斯分布和混合威布爾分布一般可以在基函數(shù)個數(shù)[1，10]和[1，5]之間找到最優(yōu)解，若在初設區(qū)間內(nèi)沒有最優(yōu)基函數(shù)個數(shù)，則區(qū)間進行倍數(shù)迭代擴大，直至找到極小值位置停止。由于實際作業(yè)工況下載荷樣本數(shù)據(jù)無法直觀呈現(xiàn)基函數(shù)個數(shù)，本文采用AIC準則（Akaike Information Criterion）與BIC準則（Bayesian Information Criterion）相結(jié)合的方法確定最優(yōu)基函數(shù)個數(shù)。AIC準則通過綜合模型擬合效果以及模型復雜程度來確定出基函數(shù)個數(shù)，對小樣本約束性更強；BIC準則在AIC基礎上引入后驗概率檢驗，對大樣本約束性更強，因此，AIC與BIC相結(jié)合的方法可先忽略樣本容量，若AIC與BIC得到的基函數(shù)個數(shù)一致，則基函數(shù)最優(yōu)個數(shù)確定，若不一致，則需要根據(jù)載荷樣本容量進一步確定基函數(shù)個數(shù)。根據(jù)AIC和BIC求解公式如式（3）～（4），求取結(jié)果中最小值所對應的值即為最優(yōu)基函數(shù)個數(shù)。

式中表示混合分布中基函數(shù)個數(shù)，表示混合分布極大對數(shù)似然函數(shù)值，為樣本數(shù)量。

設混合高斯分布基函數(shù)最優(yōu)個數(shù)取值[1，10]，混合威布爾分布基函數(shù)最優(yōu)個數(shù)取值[1，5]，對式（3）～（4）計算得到混合分布基函數(shù)個數(shù)與AIC、BIC值的曲線圖如圖4所示。

圖4 基函數(shù)個數(shù)與基于AIC、BIC準則的值 Fig.4 Number of basis functions and the values based on the AIC（Akaike Information Criterion）,BIC（Bayesian Information Criterion）

由圖4可知，當混合高斯分布基函數(shù)個數(shù)為3時，AIC與BIC值均最小，可認為均值服從基函數(shù)個數(shù)為3的混合高斯分布；當混合威布爾分布基函數(shù)個數(shù)為2時，AIC與BIC值均最小，可認為幅值服從基函數(shù)個數(shù)為2的混合威布爾分布。在實際的應用過程中，AIC準則與BIC準則得到最優(yōu)基函數(shù)個數(shù)的值可能不同，此時需要根據(jù)樣本容量的大小進一步確定基函數(shù)最優(yōu)個數(shù)。

2.3 EM算法最優(yōu)化參數(shù)估計

傳統(tǒng)參數(shù)估計方法，必須在載荷數(shù)據(jù)比較完整的情況下，才能準確估計載荷分布參數(shù)。首先載荷根據(jù)基函數(shù)個數(shù)分成組，其次將分組后載荷進行參數(shù)估計，在參數(shù)估計過程中存在最優(yōu)化分組以及部分數(shù)據(jù)未知無法直接得到分布參數(shù)等問題，因此，本文采用EM算法對載荷進行最優(yōu)化分組與參數(shù)估計。以混合高斯分布最優(yōu)化分組參數(shù)估計過程為例，過程分為E步和M步。首先，簡化混合高斯模型，并求對數(shù)似然函數(shù)，進行E步，根據(jù)當前估計參數(shù)計算屬于該類估計分布的后驗概率，進行M步，求解估計參數(shù)，通過不斷迭代計算，設第步迭代估計參數(shù)為 Θ，當參數(shù) Θ使得對數(shù)似然函數(shù)最大時，即 ‖ Θ-Θ‖足夠小時，參數(shù) Θ為最優(yōu)估計參數(shù)。

經(jīng)過EM算法對混合分布見式（2）進行優(yōu)化分組和參數(shù)估計，得到混合高斯分布與混合威布爾分布參數(shù)估計求解結(jié)果如表1所示。

表1 混合分布參數(shù)估計求解結(jié)果 Table 1 Mixed distribution parameter estimation results

2.4 混合分布模型擬合效果檢驗

為驗證擬合效果，分別從數(shù)據(jù)變化、誤差、曲線擬合三個方面進行擬合效果檢驗。結(jié)果如表2及圖5所示，從擬合效果檢驗參數(shù)看，單分布與混合分布擬合均通過檢驗，但是混合分布擬合效果明顯提高。從數(shù)據(jù)變化角度，與單分布相比，混合高斯分布決定系數(shù)、相關系數(shù)分別提高3.45%和1.73%，混合威布爾決定系數(shù)、相關系數(shù)分別提高6.02%和4.87%；在誤差角度，混合高斯分布均方根誤差與殘差平方和分別降低31.38%和52.92%，混合威布爾均方根誤差與殘差平方和分別降低40%和64.01%，檢驗結(jié)果表明混合分布更能準確描述載荷分布情況。

表2 混合高斯分布及混合威布爾分布擬合參數(shù)比較檢驗 Table 2 Mixed Gaussian and mixed Weibull fitting parameter comparison test

從概率密度擬合曲線比較圖看，混合分布概率密度曲線走勢基本與載荷數(shù)據(jù)一致，單分布在圖5a、5c標注區(qū)域無法完全覆蓋載荷頻次信息，容易出現(xiàn)欠擬合；從累積分布函數(shù)（Cumulative Distribution Function，CDF）擬合曲線比較圖看，可以認為混合分布CDF完全與實際CDF曲線吻合，而單分布與實際CDF曲線在圖5b、5d標注區(qū)域有較大偏離，表明單分布在描述實測數(shù)據(jù)的載荷分布還有較大誤差。

圖5 混合分布與單分布擬合比較 Fig.5 Comparison of mixed distribution and single distribution fitting

3 程序載荷譜編制與臺架動態(tài)加載試驗驗證

3.1 程序載荷譜編制及驗證

根據(jù)圖3a的64級雨流矩陣聯(lián)合概率密度函數(shù)可得到各外推倍數(shù)下每個均幅值區(qū)間的頻次，目前八級載荷譜最為常用，均值、幅值均按照極值與比例系數(shù)（1、0.95、0.85、0.725、0.575、0.425、0.275、0.125）的乘積進行八級分級，并根據(jù)混合高斯與混合威布爾聯(lián)合概率分布函數(shù)計算均幅值每個區(qū)間內(nèi)對應出現(xiàn)的載荷頻次，得到牽引力八級二維程序載荷譜如表3所示。

表3 牽引力八級二維程序載荷譜 Table 3 Eight-stage two-dimensional program load spectrum of traction

由于八級二維程序載荷譜難以用于加載系統(tǒng)進行加載，所以通常將二維載荷譜按照變均值法編制成一維程序載荷譜，求出均值波動中心，然后對波動中心進行施加相應級別載荷。波動中心為各分級載荷與頻次乘積之和/總頻次，得到波動中心為1 645.15N。基于Miner疲勞累積損傷理論以及S-N曲線，得到轉(zhuǎn)換前后損傷累積公式如式（5）～（6），由轉(zhuǎn)換前后疲勞損傷一致即=FD得到一維程序載荷譜的加載頻次如表4所示，根據(jù)加載頻次編制“低-高-低”的八級一維程序載荷譜，用于臺架試驗加載。

表4 牽引力八級一維程序載荷譜 Table 4 Eight-stage one-dimensional program load spectrum of traction

式中FD為二維程序載荷譜疲勞損傷；FD為一維程序載荷譜疲勞損傷；F為二維程序載荷譜分級幅值；n為的對應頻次；N為應力下疲勞壽命的總頻次；F為一維程序載荷譜分級目標幅值；n為的對應頻次；為材料性質(zhì)常數(shù)；為95%可靠度下材料對應的常數(shù)。

為更好分析外推前后的載荷數(shù)據(jù)統(tǒng)計特性以及均幅值頻次信息，外推前后數(shù)據(jù)統(tǒng)計特性如表5所示，外推一倍載荷與原始載荷均幅值頻次圖如圖6所示。從表5可以看出極值范圍變大，擴大為原來的1.194倍，均值提高，但標準差和方差基本不變；從載荷頻次圖來看，外推一倍載荷譜與原始載荷均值幅值頻次基本保持一致，均值頻次整體向右偏移對應頻次一致，幅值頻次范圍沒有變化，說明外推載荷整體變大，而波動情況基本沒有發(fā)生變化。

圖6 外推一倍載荷譜與原始載荷頻次對比 Fig.6 Comparison of extrapolated double load spectrum and original load frequency

表5 外推前后載荷統(tǒng)計特性 Table5 Statistical characteristics of data before and after extrapolation

3.2 基于PID的拖拉機三點懸掛外推載荷譜臺架試驗驗證

如圖7所示，拖拉機三點懸掛牽引力室內(nèi)臺架加載系統(tǒng)由加載油缸、懸掛提升框架、牽引力傳感器、拉繩傳感器、角度傳感器、電液比例溢流閥、比例換向閥及液壓泵站等組成，設備參數(shù)如表6所示。拖拉機三點懸掛提升框架通過牽引力傳感器與加載油缸相連，角度傳感器安裝在加載油缸支架底部；加載油缸對提升框架施加阻力，用于模擬拖拉機田間作業(yè)承受的懸掛載荷，并通過牽引力傳感器與角度傳感器檢測模擬載荷大小及加載角度。加載力、加載方向及伸縮速度通過液壓泵站、電液比例溢流閥及電液比例換向閥進行控制。此外，基于NI FPGA設計開發(fā)了拖拉機三點懸掛加載平臺測控儀，用于傳感器信號的實時采集與輸出控制，實現(xiàn)數(shù)據(jù)的顯示、保存及預處理等。

圖7 拖拉機三點懸掛牽引力加載系統(tǒng)總體方案與實物圖 Fig.7 Overall scheme and diagram of tractor three-point suspension traction loading system

表6 加載系統(tǒng)設備參數(shù) Table 6 Device parameters of loading system

三點懸掛牽引力試驗中，為簡化控制系統(tǒng)的復雜性，保持系統(tǒng)壓力與流量穩(wěn)定，通過調(diào)節(jié)電液比例溢流閥回路壓力，實現(xiàn)對加載油缸加載力的控制。試驗中，將牽引力外推載荷譜作為目標信號用于加載系統(tǒng)，牽引力傳感器實時采集加載力，基于PID反饋調(diào)節(jié)實現(xiàn)牽引力載荷譜的動態(tài)加載與復現(xiàn)。

由于加載設備動態(tài)響應性能有限，目前機械零部件疲勞壽命試驗方法通常采用靜態(tài)載荷加載或靜態(tài)逐級加載的方式，無法真實還原田間作業(yè)載荷特征，影響機械零部件產(chǎn)品設計及疲勞壽命預測精度等。本文將外推一倍載荷譜用于加載試驗，加載頻率為20 Hz。牽引力載荷譜的加載效果如圖8所示。

通過對實測牽引力和外推載荷譜輸入牽引力曲線進行分析，兩者均方差為280 N，系統(tǒng)最大超調(diào)量為5.19%，最大穩(wěn)定時間為0.2 s，最大延時為0.4 s，系統(tǒng)動態(tài)響應效果較好；從圖8可知加載系統(tǒng)存在0～450 N范圍的加載死區(qū)，由于控制系統(tǒng)加載時，電液比例換向閥方向和開口大小確定，部分液壓油經(jīng)過換向閥進入油缸所致。通過學者的研究發(fā)現(xiàn)，刪除最大載荷的10%～15%不會對零件的疲勞壽命產(chǎn)生明顯影響，外推數(shù)據(jù)的最大值為4.657 kN，刪除最大值的10%以下小載荷后，外推數(shù)據(jù)最小值為465.7 N，大于死區(qū)的范圍，因此，可以認為加載死區(qū)對外推載荷譜加載影響可以忽略。從放大曲線可以看出，牽引力基本可以實時跟隨載荷譜輸入，室內(nèi)臺架系統(tǒng)能夠滿足拖拉機三點懸掛牽引力的動態(tài)加載和復現(xiàn)。

圖8 牽引加載力與反饋力跟隨曲線 Fig.8 Loading force and feedback force curves under different loading frequency

由Miner疲勞累積損傷理論以及S-N曲線可知，經(jīng)過雨流計數(shù)之后，不同應力的次數(shù)與此應力下疲勞壽命的比值和代表該段應力數(shù)據(jù)的總壽命。對外推載荷、加載載荷進行雨流計數(shù)分析，得到均值、幅值頻次對比如圖9所示，可以看出外推一倍載荷與加載反饋載荷譜均值頻次、幅值頻次基本一致，可認為原始載荷、外推載荷、加載載荷造成的損傷等效，為基于載荷譜編制的拖拉機室內(nèi)臺架試驗提供理論依據(jù)和參考。

圖9 外推載荷譜與加載載荷譜頻次對比 Fig.9 Frequency comparison between extrapolated load spectrum and loaded load spectrum

4 結(jié) 論

1）以拖拉機三點懸掛牽引力載荷為例，建立混合分布模型，采用AIC準則與BIC準則相結(jié)合的方法確定混合分布最優(yōu)基函數(shù)個數(shù)，并基于EM算法進行最優(yōu)分布擬合求解，通過與傳統(tǒng)參數(shù)外推單分布擬合相比，混合高斯分布與混合威布爾決定系數(shù)分別提高3.45%和6.02%，相關系數(shù)分別提高1.73%和4.87%；均方根誤差分別降低31.38%和40%，殘差平方和分別降低52.92%和64.01%，驗證結(jié)果表明該方法的擬合效果總體更能準確描述載荷分布情況。

2）提出了基于最優(yōu)分布擬合的載荷譜編制方法，將原始載荷極值和頻次外推至10編制一維程序載荷譜，并將外推一倍載荷譜與原始載荷數(shù)據(jù)進行對比分析：外推后載荷范圍擴大為原來的1.194倍，標準差和方差基本保持不變；雨流計數(shù)后均值頻次整體向右偏移但對應頻次一致，幅值頻次范圍沒有變化，說明外推載荷極值范圍得到擴充，疲勞損傷進一步擴大，能夠在保留原始載荷規(guī)律特性的前提下實現(xiàn)均幅值的雙向外推。

3）開發(fā)了一套拖拉機三點懸掛牽引力室內(nèi)臺架加載系統(tǒng)，并基于PID控制實現(xiàn)懸掛牽引力的反饋調(diào)節(jié)。將基于上述方法的外推牽引力載荷譜數(shù)據(jù)與實測牽引力數(shù)據(jù)進行對比分析，結(jié)果表明：加載系統(tǒng)最大超調(diào)量為5.19%，最大穩(wěn)定時間為0.2 s，最大延時為0.4 s，室內(nèi)臺架系統(tǒng)能夠滿足拖拉機三點懸掛牽引力的加載與復現(xiàn)，并對加載載荷譜與反饋載荷譜進行雨流計數(shù)分析，均值、幅值頻次一致，兩者疲勞損傷等效，為基于載荷譜編制的拖拉機室內(nèi)臺架試驗提供理論依據(jù)和參考。