利用優(yōu)選狀態(tài)數(shù)的MCMC模擬農(nóng)機(jī)裝備負(fù)載

楊子涵，宋正河

利用優(yōu)選狀態(tài)數(shù)的MCMC模擬農(nóng)機(jī)裝備負(fù)載

楊子涵，宋正河※

（中國(guó)農(nóng)業(yè)大學(xué)現(xiàn)代農(nóng)業(yè)裝備優(yōu)化設(shè)計(jì)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100083）

傳統(tǒng)馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carlo，MCMC）方法中狀態(tài)數(shù)的選取常依賴于主觀經(jīng)驗(yàn)，用于農(nóng)機(jī)裝備負(fù)載模擬時(shí)，狀態(tài)數(shù)取值不當(dāng)將導(dǎo)致負(fù)載模擬精度降低或算法運(yùn)行時(shí)間冗長(zhǎng)。針對(duì)此問(wèn)題，該研究提出一種基于偽損傷一致性的狀態(tài)數(shù)優(yōu)選方法。首先確定MCMC算法中狀態(tài)數(shù)的初選范圍，然后分別計(jì)算范圍內(nèi)不同狀態(tài)數(shù)所對(duì)應(yīng)的負(fù)載模擬結(jié)果，最后以生成的模擬負(fù)載與原始載荷之間的損傷一致性為評(píng)價(jià)準(zhǔn)則確定優(yōu)選狀態(tài)數(shù)。利用拖拉機(jī)關(guān)鍵零部件的實(shí)測(cè)載荷數(shù)據(jù)對(duì)該方法進(jìn)行驗(yàn)證。結(jié)果表明，隨著狀態(tài)數(shù)的提高，模擬負(fù)載與原始載荷之間的損傷一致性變化趨于平穩(wěn)，算法運(yùn)算時(shí)長(zhǎng)增速不斷提高，相比于傳統(tǒng)方法，基于優(yōu)選狀態(tài)數(shù)的MCMC算法能夠得到偽損傷差異在1%以內(nèi)的負(fù)載模擬結(jié)果，與載荷譜編制的目標(biāo)需求更加匹配，在保證模擬結(jié)果精度的同時(shí)有效減少運(yùn)算成本。該研究能夠?yàn)檗r(nóng)機(jī)裝備關(guān)鍵零部件的動(dòng)態(tài)仿真分析及可靠性試驗(yàn)提供更加可靠的數(shù)據(jù)支撐。

農(nóng)業(yè)機(jī)械；模擬；載荷；馬爾科夫鏈；蒙特卡洛法；優(yōu)選狀態(tài)數(shù)；偽損傷

0 引 言

隨著農(nóng)業(yè)裝備不斷向現(xiàn)代化、智能化和規(guī)?；l(fā)展，農(nóng)業(yè)裝備生產(chǎn)企業(yè)愈發(fā)重視產(chǎn)品的可靠性問(wèn)題[1-2]。如今，為了爭(zhēng)奪更廣闊的市場(chǎng)資源，企業(yè)必須以適應(yīng)不同地域、不同作業(yè)需求為目標(biāo)生產(chǎn)耐用可靠的機(jī)器[3]。為實(shí)現(xiàn)這一目標(biāo)，研發(fā)人員需在產(chǎn)品投放市場(chǎng)前能夠?qū)Ξa(chǎn)品整機(jī)及關(guān)鍵零部件的疲勞壽命做出更為精準(zhǔn)的預(yù)估[4-5]。動(dòng)力機(jī)械的載荷不穩(wěn)定、隨機(jī)波動(dòng)和隨時(shí)間變化等特征是對(duì)其可靠性產(chǎn)生威脅的重要因素[6]。因此，如何準(zhǔn)確獲取產(chǎn)品在實(shí)際作業(yè)工況下的載荷特征，編制全生命周期載荷譜成為了可靠性研究領(lǐng)域的一個(gè)重要課題[7]。相比于一般動(dòng)力機(jī)械，農(nóng)業(yè)裝備的作業(yè)環(huán)境更加惡劣，工況復(fù)雜多變且具有明顯的季節(jié)性作業(yè)特征，使農(nóng)業(yè)裝備載荷譜編制變得更加復(fù)雜。

針對(duì)飛機(jī)關(guān)鍵部件結(jié)構(gòu)試驗(yàn)的載荷譜最先被提出，這種譜塊形式的程序載荷譜主要以疲勞線性累計(jì)理論為支撐，對(duì)于臺(tái)架動(dòng)態(tài)加載性能要求較低，因此被廣泛應(yīng)用于車輛、船舶、風(fēng)機(jī)等領(lǐng)域[8-9]。隨著研究的深入，越來(lái)越多的文獻(xiàn)表明，載荷的均值、加載的順序和頻率同樣會(huì)對(duì)結(jié)構(gòu)的疲勞失效產(chǎn)生不可忽略影響[10]。因此，更為細(xì)膩且考慮載荷加載順序的載荷譜編制方法開(kāi)始出現(xiàn)，如考慮變均值混合工況的車輛載荷譜CARLOS和海洋平臺(tái)載荷譜Wash[11]。國(guó)內(nèi)閆楚良等[12]根據(jù)飛機(jī)實(shí)際受載特征提出了高置信度中值隨機(jī)疲勞載荷譜編制方法。隨著電氣和液壓伺服系統(tǒng)的出現(xiàn)及有限元仿真技術(shù)的快速發(fā)展，如今，載荷譜的加載形式已不再局限于簡(jiǎn)化后的譜塊形式[13]，在此背景下，如何基于短時(shí)的實(shí)測(cè)載荷信息模擬出能夠滿足試驗(yàn)臺(tái)長(zhǎng)期加載及全生命周期仿真需求的載荷時(shí)間序列，進(jìn)一步縮減測(cè)試成本成為了新的研究熱點(diǎn)。

目前已有的負(fù)載模擬方法主要有雨流域、時(shí)域、頻域等。雨流域的相關(guān)研究[14-15]發(fā)展較早且技術(shù)最為成熟，但是由于無(wú)法保留負(fù)載的時(shí)序信息限制了其應(yīng)用范圍。從頻域角度對(duì)負(fù)載進(jìn)行建模能夠反映出實(shí)測(cè)載荷樣本的物理含義，目前應(yīng)用對(duì)象主要集中于頻譜特征明顯的旋轉(zhuǎn)部件[16-17]。近年來(lái)又有學(xué)者將機(jī)器學(xué)習(xí)的相關(guān)方法應(yīng)用于載荷模擬[18]，由于對(duì)載荷樣本要求較高并且需人工劃分工況，因此未得到廣泛應(yīng)用。

以獲取與實(shí)測(cè)信號(hào)時(shí)域特征相一致的載荷時(shí)間序列為目標(biāo)的時(shí)域模擬技術(shù)是目前的主流研究方向。Johannesson[19]提出了基于極值理論載荷時(shí)域外推方法，隨后眾多學(xué)者針對(duì)該方法的關(guān)鍵技術(shù)進(jìn)行了深入研究，然而該類方法僅對(duì)超出閾值的極端載荷進(jìn)行建模，忽略了對(duì)一般載荷的模擬重構(gòu)。MCMC（Markov chain Monte Carlo）方法能夠基于負(fù)載變化特征更全面地模擬載荷時(shí)間序列。Rychlik等[20]提出了利用馬爾科夫鏈對(duì)隨機(jī)負(fù)載進(jìn)行建模的方法，并利用實(shí)測(cè)卡車負(fù)載數(shù)據(jù)對(duì)該方法進(jìn)行了驗(yàn)證。Carboni等[21]利用馬爾科夫狀態(tài)轉(zhuǎn)移矩陣對(duì)液壓泵負(fù)載進(jìn)行建模，并利用負(fù)載累積時(shí)間曲線重構(gòu)負(fù)載的持續(xù)時(shí)長(zhǎng)，解決了馬爾科夫域無(wú)法保留時(shí)域信息的問(wèn)題。武玉倩[22]討論了駕駛員行為特征對(duì)負(fù)載的影響，提出了基于隱馬爾科夫模型的載荷譜外推方法。Cerrini等[23]針對(duì)混凝土灌裝機(jī)械的鉸接臂負(fù)載，提出了一種改進(jìn)的馬爾科夫建模方法，能夠考慮狀態(tài)切換所造成的負(fù)載差異。Wang等[24]基于MCMC方法提出了應(yīng)用于裝載機(jī)非平穩(wěn)負(fù)載的循環(huán)工況模擬方法。劉純等[25]將MCMC方法應(yīng)用于風(fēng)電出力時(shí)間序列的模擬。在農(nóng)機(jī)裝備領(lǐng)域，Paraforos等[26-27]利用馬爾科夫鏈對(duì)摟草機(jī)械機(jī)體應(yīng)變載荷和拖拉機(jī)油耗進(jìn)行建模，通過(guò)蒙特卡洛方法生成了時(shí)頻特征較為一致的模擬結(jié)果。

需要注意的是，上述MCMC方法中馬爾科夫鏈的狀態(tài)數(shù)選取均依賴于主觀經(jīng)驗(yàn)，如文獻(xiàn)[24]中所描述的狀態(tài)數(shù)取值一般為32或64，文獻(xiàn)[27]中選取的狀態(tài)數(shù)則為256。徐沈智等[28]以模擬序列與原始序列的自相關(guān)特性為評(píng)判準(zhǔn)則，討論了風(fēng)電功率序列建模時(shí)狀態(tài)數(shù)對(duì)于模擬結(jié)果的影響，并提出了基于ACF（Auto-Correlation Function）曲線距離最小的狀態(tài)數(shù)選取方法。然而風(fēng)電功率序列與農(nóng)機(jī)裝備負(fù)載在數(shù)據(jù)特征與應(yīng)用方式上有著本質(zhì)差異，因此并不適用農(nóng)機(jī)裝備負(fù)載模擬。

基于此，本文首先討論了MCMC方法應(yīng)用于農(nóng)機(jī)裝備負(fù)載模擬時(shí)，狀態(tài)數(shù)取值對(duì)于模擬結(jié)果的影響。在此基礎(chǔ)上，以保障載荷譜編制精度為目標(biāo)，提出了一種基于偽損傷一致性的狀態(tài)數(shù)優(yōu)選方法。最后以拖拉機(jī)犁耕工況下前橋振動(dòng)負(fù)載為例驗(yàn)證了狀態(tài)數(shù)選取的合理性，進(jìn)一步討論了農(nóng)機(jī)裝備不同作業(yè)段負(fù)載差異對(duì)于狀態(tài)數(shù)選取結(jié)果的影響，并給出了基于優(yōu)選狀態(tài)數(shù)的MCMC算法在農(nóng)機(jī)裝備負(fù)載模擬中的應(yīng)用范例。

1 基于優(yōu)選狀態(tài)數(shù)的MCMC方法

1.1 MCMC方法

需要注意的是，與傳統(tǒng)馬爾科夫鏈應(yīng)用方式不同，在對(duì)實(shí)測(cè)載荷序列進(jìn)行建模時(shí)，往往只關(guān)注與局部極端載荷，即實(shí)測(cè)載荷序列中局部最大值與最小值所組成的轉(zhuǎn)折點(diǎn)序列，這種方式能夠在保留載荷時(shí)域信息的同時(shí)有效縮減數(shù)據(jù)量。利用MCMC方法對(duì)實(shí)測(cè)載荷進(jìn)行模擬的具體步驟如下：

1）載荷時(shí)域特征提取。首先對(duì)實(shí)測(cè)載荷信號(hào)進(jìn)行峰谷值提取，得到轉(zhuǎn)折點(diǎn)序列{X}，然后將峰值序列{1,2,…}與谷值序列{1,2,…}分別進(jìn)行提取，并記錄升程S對(duì)應(yīng)的時(shí)間間隔ΔT和降程s對(duì)應(yīng)的時(shí)間間隔Δt：

為便于蒙特卡洛模擬，進(jìn)一步計(jì)算上述概率矩陣的累計(jì)概率矩陣：

1.2 狀態(tài)數(shù)選取對(duì)載荷模擬結(jié)果的影響

本文選取拖拉機(jī)犁耕作業(yè)和整地作業(yè)工況下關(guān)鍵零部件的振動(dòng)、應(yīng)力和扭矩負(fù)載作為驗(yàn)證數(shù)據(jù)集，如圖1所示。

犁耕工況下前橋的振動(dòng)和應(yīng)力載荷數(shù)據(jù)采集于北京上莊試驗(yàn)站，拖拉機(jī)型號(hào)為約翰迪爾6B-1354，掛接機(jī)具為1LF-535型液壓翻轉(zhuǎn)犁，犁耕作業(yè)時(shí)拖拉機(jī)檔位為B2擋，作業(yè)速度為5～10 km/h，作業(yè)路徑如圖2a所示。整地工況下的傳動(dòng)軸扭矩載荷數(shù)據(jù)采集于江蘇鹽城新洋試驗(yàn)站，拖拉機(jī)型號(hào)為東方紅LX2204，掛接機(jī)具為HR4004D型驅(qū)動(dòng)靶，整地作業(yè)時(shí)拖拉機(jī)檔位為中三檔，作業(yè)速度為2～6 km/h，作業(yè)路徑如圖2b所示。

上述兩種工況下的田間試驗(yàn)均于當(dāng)?shù)卮焊鳂I(yè)季內(nèi)開(kāi)展，試驗(yàn)嚴(yán)格參照標(biāo)準(zhǔn)T/NJ 1252—2020《拖拉機(jī)田間作業(yè)機(jī)組試驗(yàn)規(guī)程犁耕作業(yè)》和T/NJ 1253—2020《拖拉機(jī)田間作業(yè)機(jī)組試驗(yàn)規(guī)程整地作業(yè)》進(jìn)行，試驗(yàn)用傳感器及設(shè)備信息如表1所示。選用當(dāng)?shù)亟?jīng)驗(yàn)豐富的機(jī)手駕駛拖拉機(jī)，試驗(yàn)過(guò)程中不存在對(duì)于機(jī)手作業(yè)方式的干預(yù)，采集的載荷數(shù)據(jù)能夠有效代表當(dāng)?shù)貙?shí)際的農(nóng)業(yè)生產(chǎn)作業(yè)特點(diǎn)。

表1 傳感器及設(shè)備信息

由圖1和圖2可知，田間作業(yè)時(shí)隨著拖拉機(jī)行駛狀態(tài)的循環(huán)切換，作業(yè)機(jī)具與拖拉機(jī)之間、作業(yè)機(jī)具與土壤之間的力和位置關(guān)系發(fā)生改變，拖拉機(jī)關(guān)鍵零部件所承受的負(fù)載特征也同步產(chǎn)生變化。以犁耕工況為例，根據(jù)拖拉機(jī)行駛狀態(tài)，可以將該工況細(xì)分為作業(yè)階段和調(diào)整階段，兩個(gè)階段交替切換構(gòu)成目標(biāo)地塊的一次完整作業(yè)，如圖3a所示。

圖3b中根據(jù)拖拉機(jī)行駛狀態(tài)不同將實(shí)測(cè)載荷劃分為16個(gè)載荷片段，可以看出，拖拉機(jī)犁耕工況負(fù)載在作業(yè)階段（O1～O8）和調(diào)整階段（A1～A8）具有明顯的特征差異。作業(yè)階段拖拉機(jī)牽引負(fù)荷較大，當(dāng)犁具調(diào)整完畢后耕作深度保持相對(duì)穩(wěn)定，因此振動(dòng)量級(jí)偏大但是各負(fù)載之間差異較小。調(diào)整階段機(jī)具處于提升狀態(tài)，牽引阻力較小，振動(dòng)主要由路面激勵(lì)引起，相對(duì)于作業(yè)階段載荷量級(jí)較小，但是由于田間路面及車速差異明顯，波動(dòng)較為劇烈且負(fù)載之間差異明顯。

MCMC方法在處理由特征差異明顯的片段所組成的載荷時(shí)間序列時(shí)，會(huì)破壞原有特征片段的獨(dú)立性。簡(jiǎn)單的解決方式即分別提取原始負(fù)載中各個(gè)特征片段分別進(jìn)行模擬，再通過(guò)循環(huán)仿真的方式實(shí)現(xiàn)負(fù)載重構(gòu)[24]。因此，本文首先以拖拉機(jī)犁耕工況下前橋橋臂所受的振動(dòng)載荷為研究對(duì)象，截取單次作業(yè)片段的載荷數(shù)據(jù)進(jìn)行分析，討論狀態(tài)數(shù)選取對(duì)載荷模擬結(jié)果的影響。

可靠性常用雨流計(jì)數(shù)法得到的雨流矩陣來(lái)統(tǒng)計(jì)載荷循環(huán)的特征信息，因此可用雨流矩陣間的相似程度來(lái)間接評(píng)價(jià)模擬結(jié)果與實(shí)測(cè)載荷的一致性。雨流差值矩陣的-范數(shù)（Frobenius norm）可由下式得出：

由圖4a可知，隨著狀態(tài)數(shù)增加，MCMC方法模擬結(jié)果的均值相對(duì)誤差和標(biāo)準(zhǔn)差相對(duì)誤差會(huì)迅速降低并趨于穩(wěn)定，相比于均值相對(duì)誤差，標(biāo)準(zhǔn)差相對(duì)誤差會(huì)更快趨于穩(wěn)定。與宏觀的統(tǒng)計(jì)特征不同，-范數(shù)能夠從載荷循環(huán)均幅值所構(gòu)成的二維空間對(duì)模擬結(jié)果進(jìn)行相似性評(píng)價(jià)，評(píng)價(jià)結(jié)果也更為細(xì)致。由圖4b可知，隨狀態(tài)數(shù)增加，-范數(shù)逐漸降低，雖然數(shù)值下降的速度逐漸平緩但具有持續(xù)下降的趨勢(shì)。

綜上可以看出，狀態(tài)數(shù)的增加會(huì)提高M(jìn)CMC方法負(fù)載模擬的精度，但是上述特征評(píng)價(jià)指標(biāo)均不能給出客觀一致且具有實(shí)際理論意義的狀態(tài)數(shù)選取方案?；诖耍疚奶岢鲆环N基于偽損傷一致性的狀態(tài)數(shù)優(yōu)選方法。

1.3 基于偽損傷一致性的狀態(tài)數(shù)選取方法

針對(duì)農(nóng)機(jī)裝備關(guān)鍵零部件實(shí)測(cè)載荷進(jìn)行模擬的目的是指導(dǎo)產(chǎn)品的可靠性設(shè)計(jì)，因此本文嘗試從損傷角度入手指導(dǎo)MCMC方法應(yīng)用于負(fù)載模擬時(shí)狀態(tài)數(shù)的選取。構(gòu)建農(nóng)機(jī)裝備關(guān)鍵零部件載荷譜所需要考慮的載荷種類不僅僅局限于應(yīng)變或應(yīng)力載荷，還包括扭矩、壓力、振動(dòng)、牽引力等。此外，大多數(shù)農(nóng)機(jī)裝備關(guān)鍵零部件并不具備精準(zhǔn)的疲勞應(yīng)力循環(huán)曲線（S-N曲線）信息，精準(zhǔn)的疲勞壽命預(yù)測(cè)往往需要大量的試驗(yàn)，耗時(shí)費(fèi)力。

偽損傷理論能夠在一定程度上反映外部激勵(lì)載荷對(duì)結(jié)構(gòu)件造成損傷的潛在能力已被證實(shí)[29]。通過(guò)計(jì)算偽損傷，能夠在零部件S-N曲線未知的情況下，實(shí)現(xiàn)不同載荷片段之間損傷程度的定量評(píng)價(jià)。因此，本文引入偽損傷理論，通過(guò)計(jì)算載荷模擬前后的偽損傷值實(shí)現(xiàn)狀態(tài)數(shù)的定量選取?；趥螕p傷一致性的狀態(tài)數(shù)選取方法過(guò)程如下：

1）利用雨流計(jì)數(shù)法提取實(shí)測(cè)載荷中每個(gè)負(fù)載循環(huán)信息，基于修正的Miner準(zhǔn)則，將所有載荷循環(huán)的損傷值進(jìn)行累加：

式中為總累計(jì)損傷；1/δ為實(shí)測(cè)第個(gè)載荷循環(huán)所造成的損傷值；為載荷循環(huán)次數(shù)；為結(jié)構(gòu)材料S-N曲線的反斜率系數(shù)；S為第個(gè)載荷循環(huán)幅值；為常數(shù)變量。

2）由于偽損傷與結(jié)構(gòu)材料參數(shù)并無(wú)嚴(yán)格對(duì)應(yīng)關(guān)系，因此偽損傷簡(jiǎn)化為

LDH檢測(cè)：空腹采集患者靜脈血，不抗凝，離心分離血清。應(yīng)用LDH檢測(cè)試劑盒（成都元和華盛股份有限公司產(chǎn)品）和i2000全自動(dòng)生化分析儀（美國(guó)Abbott Laboratories公司產(chǎn)品）檢測(cè)LDH。正常值參考范圍：150～245 U/L。

由公式（12）可知，偽損傷僅與載荷循環(huán)幅值和反映材料疲勞特性的反斜率系數(shù)有關(guān)。借鑒汽車零部件的偽損傷計(jì)算準(zhǔn)則，對(duì)于焊接結(jié)構(gòu)件或發(fā)生裂紋擴(kuò)展的零部件，取3；對(duì)于典型的表面粗糙零部件，取5；對(duì)于具有光滑表面的零部件，取7[29]。

3）分別計(jì)算不同狀態(tài)數(shù)下模擬負(fù)載產(chǎn)生的偽損傷，進(jìn)一步利用公式（13）計(jì)算不同狀態(tài)數(shù)對(duì)應(yīng)的偽損傷系數(shù)。越接近于1，模擬負(fù)載與原始負(fù)載的偽損傷一致性越高。

式中0為原始負(fù)載對(duì)應(yīng)的偽損傷。

設(shè)定偽損傷系數(shù)的上限閾值為Q，選取低于閾值且數(shù)值最小的狀態(tài)數(shù)作為優(yōu)選狀態(tài)數(shù)。綜上，利用優(yōu)選狀態(tài)數(shù)的MCMC方法流程如圖5所示。

2 結(jié)果與分析

為驗(yàn)證本文狀態(tài)數(shù)選取方法的有效性，選取1.2節(jié)中截取的單程犁耕作業(yè)片段進(jìn)行實(shí)例分析。首先，依次計(jì)算不同狀態(tài)數(shù)取值所對(duì)應(yīng)的原始載荷等長(zhǎng)度模擬結(jié)果的偽損傷系數(shù)。本文狀態(tài)數(shù)的初選范圍為[4,300]，由于負(fù)載模擬過(guò)程具有一定隨機(jī)性，每次模擬結(jié)果的偽損傷系數(shù)并不完全一致。為消除隨機(jī)性對(duì)于負(fù)載模擬結(jié)果的影響，針對(duì)每個(gè)狀態(tài)數(shù)取值分別進(jìn)行100次獨(dú)立重復(fù)試驗(yàn)，計(jì)算結(jié)果取平均值作為該狀態(tài)數(shù)對(duì)應(yīng)偽損傷系數(shù)的最終結(jié)果。狀態(tài)數(shù)與偽損傷系數(shù)的對(duì)應(yīng)關(guān)系及算法用時(shí)如圖6所示。

由圖6可知，與原始載荷相比，MCMC方法得到的載荷模擬結(jié)果的偽損傷系數(shù)整體偏大，隨著狀態(tài)數(shù)增加，模擬結(jié)果的偽損傷系數(shù)會(huì)迅速降低后趨近與1，而算法的運(yùn)算時(shí)長(zhǎng)保持有上升的趨勢(shì)且增速不斷提高。為在滿足模擬精度的前提下盡可能縮減運(yùn)算成本，本文將偽損傷系數(shù)的閾值設(shè)定為1.01，即與原始載荷偽損傷的差異在1%以內(nèi)，得出該負(fù)載片段的優(yōu)選狀態(tài)數(shù)為34。

將原始載荷與模擬結(jié)果分別進(jìn)行雨流計(jì)數(shù)并繪制載荷循環(huán)幅值累計(jì)頻次曲線，如圖7所示?？梢钥闯?，狀態(tài)數(shù)取值過(guò)低時(shí)，模擬負(fù)載的循環(huán)幅值偏大，致使生成的載荷譜更為惡劣，在狀態(tài)數(shù)達(dá)到優(yōu)選狀態(tài)數(shù)后，比優(yōu)選狀態(tài)數(shù)更大的取值對(duì)提高載荷模擬精度的作用變得十分有限，這和偽損傷系數(shù)與狀態(tài)數(shù)的變化關(guān)系相一致，驗(yàn)證了本文狀態(tài)數(shù)選取方法的有效性。

按照本文方法對(duì)1.2節(jié)各載荷片段（O1～O8，A1～A8）依次進(jìn)行分析，計(jì)算出各片段的優(yōu)選狀態(tài)數(shù)及相應(yīng)的偽損傷系數(shù)，如表2所示。由表2可知，基于偽損傷一致性的狀態(tài)數(shù)選取方法可以有效適用于犁耕工況下各振動(dòng)載荷片段，基于優(yōu)選狀態(tài)數(shù)的負(fù)載模擬結(jié)果均能夠滿足設(shè)定的偽損傷一致性要求。此外，受農(nóng)機(jī)裝備田間作業(yè)特點(diǎn)影響，即使在單一犁耕工況下，由于載荷特征不一致，不同載荷片段對(duì)應(yīng)的優(yōu)選狀態(tài)數(shù)計(jì)算結(jié)果同樣具有差異性。相比于作業(yè)階段，調(diào)整階段各載荷片段之間的差異更為明顯，優(yōu)選狀態(tài)數(shù)計(jì)算結(jié)果也更加分散。因此，在利用MCMC方法對(duì)農(nóng)機(jī)裝備作業(yè)載荷進(jìn)行模擬時(shí)有必要針對(duì)不同負(fù)載片段分別進(jìn)行優(yōu)選狀態(tài)數(shù)的確定。

表2 犁耕工況下各載荷片段的優(yōu)選狀態(tài)數(shù)

同理，對(duì)1.2節(jié)中犁耕工況下前橋應(yīng)力負(fù)載和整地工況下傳動(dòng)軸扭矩負(fù)載分別進(jìn)行循環(huán)模擬。為進(jìn)一步驗(yàn)證所述方法對(duì)于拖拉機(jī)各關(guān)鍵零部件在不同作業(yè)工況下所受負(fù)載的適用性，提取循環(huán)模擬結(jié)果中各載荷片段的均值、標(biāo)準(zhǔn)差和載荷循環(huán)幅值最大值作為統(tǒng)計(jì)特征值，按照負(fù)載類型不同進(jìn)行劃分，將各載荷片段模擬結(jié)果特征值與原始載荷相應(yīng)特征值的差值作為一個(gè)集合，得到各統(tǒng)計(jì)特征值的偏差范圍[min,max]如表3所示，min和max分別為差值集合中的最小值和最大值。

表3 各載荷片段模擬結(jié)果的統(tǒng)計(jì)特征值偏差范圍

作為一種基于統(tǒng)計(jì)學(xué)原理的建模方式，理論上本文方法對(duì)于不同工況下不同類型和特征的負(fù)載數(shù)據(jù)具有通用性。由表3各統(tǒng)計(jì)特征值的對(duì)比結(jié)果可知，各負(fù)載片段的模擬結(jié)果與原始負(fù)載間的特征值差值范圍很小，進(jìn)一步驗(yàn)證了本文方法用于以載荷譜編制為目標(biāo)的農(nóng)機(jī)裝備負(fù)載模擬的通用性。

3 結(jié) 論

針對(duì)傳統(tǒng)MCMC方法中狀態(tài)數(shù)選取缺乏理論依據(jù)的問(wèn)題，以拖拉機(jī)前轉(zhuǎn)向驅(qū)動(dòng)橋犁耕作業(yè)時(shí)所受到的振動(dòng)載荷為對(duì)象，研究不同狀態(tài)數(shù)取值對(duì)于負(fù)載模擬精度的影響并提出狀態(tài)數(shù)優(yōu)選方法。主要結(jié)論如下：

1）狀態(tài)數(shù)取值不同會(huì)對(duì)負(fù)載模擬結(jié)果精度產(chǎn)生重要影響。隨著狀態(tài)數(shù)取值的增加，模擬負(fù)載與原始載荷之間的均值誤差、標(biāo)準(zhǔn)差誤差和雨流矩陣偏差會(huì)快速減小并趨于穩(wěn)定。

2）提出了基于偽損傷一致性的狀態(tài)數(shù)選取方法。隨著狀態(tài)數(shù)提高，模擬負(fù)載與原始載荷之間的損傷一致性逐漸提高并趨于平穩(wěn)，而算法運(yùn)算時(shí)長(zhǎng)的增速依然不斷提高。通過(guò)設(shè)定偽損傷系數(shù)的閾值，即可計(jì)算出滿足損傷一致性要求且運(yùn)算時(shí)長(zhǎng)最小的優(yōu)選狀態(tài)數(shù)。

3）實(shí)例分析結(jié)果表明，基于優(yōu)選狀態(tài)數(shù)的MCMC方法能夠得到偽損傷差異在1%以內(nèi)的載荷模擬結(jié)果，各負(fù)載片段的模擬結(jié)果特征值與原始載荷特征值具有較高的一致性。與傳統(tǒng)方法相比，基于優(yōu)選狀態(tài)數(shù)的MCMC方法與載荷譜編制的目標(biāo)更加匹配，能夠滿足拖拉機(jī)各關(guān)鍵零部件在不同工況下的負(fù)載模擬需求，在保證負(fù)載模擬精度的前提下有效減少了運(yùn)算成本。

[1] 謝斌，武仲斌，毛恩榮. 農(nóng)業(yè)拖拉機(jī)關(guān)鍵技術(shù)發(fā)展現(xiàn)狀與展望[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2018，49(8)：1-17.

Xie Bin, Wu Zhongbin, Mao Enrong. Development and prospect of key technologies on agriculture tractor[J]. Transactions of the Chinese Society for Agricultural Machinery, 2015, 51(8): 131-136. (in Chinese with English abstract)

[2] 尹宜勇，付寧善，廖頻，等. 基于DTW距離的拖拉機(jī)傳動(dòng)軸載荷樣本長(zhǎng)度計(jì)算方法[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2020，36(12)：54-60.

Yin Yiyong, Fu Ningshan, Liao Pin, et al. Calculation method of load sample size for tractor drive shafts based on dynamic time warping distance[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(12): 54-60. (in Chinese with English abstract)

[3] Paraforos D S, Griepentrog H W, Vougioukas S G. Modeling and simulation of a four-rotor rake loading for predicting accumulated fatigue damage: A Markov regime-switching approach[J]. Applied Engineering in Agriculture, 2018, 34(2): 317-325.

[4] Wen C K, Xie B, Song Z H, et al. Methodology for designing tractor accelerated structure tests for an indoor drum-type test bench[J]. Biosystems Engineering, 2021, 205: 1-26.

[5] 楊子涵，宋正河，尹宜勇，等. 基于POT模型的大功率拖拉機(jī)傳動(dòng)軸載荷時(shí)域外推方法[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2019，35(15)：40-47.

Yang Z, Song Z, Yin Y, et al. Time domain extrapolation method for load of drive shaft of high-power tractor based on POT model[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(15): 40-47. (in Chinese with English abstract)

[6] 軒福貞，朱明亮，王國(guó)彪. 結(jié)構(gòu)疲勞百年研究的回顧與展望[J]. 機(jī)械工程學(xué)報(bào)，2021，57(6)：26-51.

Xuan Fuzhen, Zhu Mingliang, Wang Guobiao. Retrospect and prospect on century-long research of structural fatigue[J]. Journal of Mechanical Engineering, 2021, 57(6): 26-51. (in Chinese with English abstract)

[7] 陳道云，孫守光，李強(qiáng). 高速列車載荷譜推斷及擴(kuò)展方法研究[J]. 機(jī)械工程學(xué)報(bào)，2018，54(10)：151-155.

Chen Daoyun, Sun Shouguang, Li Qiang. Study on deduction and extend of high-speed train load spectrum[J]. Journal of Mechanical Engineering, 2018, 54(10): 151-155. (in Chinese with English abstract)

[8] 高云凱，徐成民，方劍光. 車身臺(tái)架疲勞試驗(yàn)程序載荷譜研究[J]. 機(jī)械工程學(xué)報(bào)，2014，50(4)：92-98.

Gao Yunkai, Xu Chengmin, Fang Jianguang. Study on the programed load spectrum of the body fatigue bench test[J]. Journal of Mechanical Engineering, 2014, 50(4): 92-98. (in Chinese with English abstract)

[9] 張英爽，王國(guó)強(qiáng)，王繼新，等. 工程車輛傳動(dòng)系載荷譜編制方法[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2011，27(4)：179-183.

Zhang Yingshuang, Wang Guoqiang, Wang Jixin, et al. Compilation method of power train load spectrum of engineering vehicle[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2011, 27(4): 179-183. (in Chinese with English abstract)

[10] Branigan Z, Suh C S. Power density-an alternative approach to quantifying fatigue failure[J]. Journal of Vibration Testing and System Dynamics, 2018, 2(4): 307-326.

[11] Heuler P, Kl?tschke H. Generation and use of standardised load spectra and load-time histories[J]. International Journal of Fatigue, 2005, 27(8): 974-990.

[12] 閻楚良，高鎮(zhèn)同. 飛機(jī)高置信度中值隨機(jī)疲勞載荷譜的編制原理[J]. 航空學(xué)報(bào)，2000(2)：118-123.

Yan Chuliang, Gao Zhentong. Compilation theory of median stochastic fatigue load spectrum with high confidence level for airplane[J]. Acta Aeronautica et Astronautica Sinica, 2000(2): 118-123. (in Chinese with English abstract)

[13] Sonsino C M. Fatigue testing under variable amplitude loading[J]. International Journal of Fatigue, 2007, 29(6): 1080-1089.

[14] Nagode M, Klemenc J, Fajdiga M. Parametric modelling and scatter prediction of rainflow matrices[J]. International Journal of Fatigue, 2001, 23(6): 525-532.

[15] Johannesson P, Thomas J J. Extrapolation of rainflow matrices[J]. Extremes, 2001, 4(3): 241-262.

[16] Poirier M, Gagnon M, Tahan A, et al. Extrapolation of dynamic load behaviour on hydroelectric turbine blades with cyclostationary modelling[J]. Mechanical Systems and Signal Processing, 2017, 82: 193-205.

[17] Wen Y K, Gu P. Description and simulation of nonstationary processes based on Hilbert spectra[J]. Journal of Engineering Mechanics, 2004, 130(8): 942-951.

[18] Qin C, Shi G, Tao J, et al. Precise cutterhead torque prediction for shield tunneling machines using a novel hybrid deep neural network[J]. Mechanical Systems and Signal Processing, 2021, 151: 107386.

[19] Johannesson P. Extrapolation of load histories and spectra[J]. Fatigue & Fracture of Engineering Materials & Structures, 2006, 29(3): 209-217.

[20] Rychlik I. Simulation of load sequences from rainflow matrices: Markov method[J]. International Journal of Fatigue, 1996, 18(7): 429-438.

[21] Carboni M, Cerrini A, Johannesson P, et al. Load spectra analysis and reconstruction for hydraulic pump components[J]. Fatigue & Fracture of Engineering Materials & Structures, 2008, 31(3‐4): 251-261.

[22] 武玉倩. 基于隱馬爾科夫模型的裝載機(jī)載荷譜編制方法研究[D]. 長(zhǎng)春：吉林大學(xué)，2017.

Wu Yuqian. Research on Load Spectra Compilation Method of Loaders Based on Hidden Markov Model[D]. Changchun: Jilin University, 2019. (in Chinese with English abstract)

[23] Cerrini A, Johannesson P, Beretta S. Superposition of manoeuvres and load spectra extrapolation[C]. Applied Mechanics and Materials, Trans Tech Publications Ltd, 2006, 5: 255-262.

[24] Wang J, Zhang J, Liang Y, et al. A cyclic simulation approach for the generation of the non-stationary load histories of engineering vehicles[J]. Journal of Mechanical Science and Technology, 2012, 26(5): 1547-1554.

[25] 劉純，呂振華，黃越輝，等. 長(zhǎng)時(shí)間尺度風(fēng)電出力時(shí)間序列建模新方法研究[J]. 電力系統(tǒng)保護(hù)與控制，2013，41(1)：7-13.

Liu Chun, Lü Zhenhua, Huang Yuehui, et al. A new method to simulate wind power time series of large time scale[J]. Power System Protection and Control, 2013, 41(1): 7-13. (in Chinese with English abstract)

[26] Paraforos D S, Griepentrog H W. Switching Markov chains for modelling the loads of a four-rotor swather under different operating modes[J]. IFAC-PapersOnLine, 2017, 50(1): 5392-5397.

[27] Paraforos D S, Griepentrog H W. Tractor fuel rate modeling and simulation using switching Markov chains on CAN-Bus data[J]. IFAC-PapersOnLine, 2019, 52(30): 379-384.

[28] 徐沈智，艾小猛，鄒佳芯，等. 優(yōu)選狀態(tài)數(shù)的MCMC算法在風(fēng)電功率序列生成中的應(yīng)用[J]. 電力自動(dòng)化設(shè)備，2019，39(5)：61-68.

Xu Shenzhi, Ai Xiaomeng, Zou Jiaxin, et al. Application of optimizing state number Markov chain Monte Carlo algorithm in wind power generation[J]. Electric Power Automation Equipment, 2019, 39(5): 61-68. (in Chinese with English abstract)

[29] 于佳偉，鄭松林，馮金芝，等. 某轎車前副車架服役載荷模擬試驗(yàn)加速方法研究[J]. 機(jī)械工程學(xué)報(bào)，2016，52(22)：112-120.

Yu Jiawei, Zheng Songlin, Feng Jinzhi, et al. Research on accelerated testing method for the service-simulation fatigue test of automotive front sub-frame[J]. Journal of Mechanical Engineering, 2016, 52(22): 112-120. (in Chinese with English abstract)

Simulation of agricultural equipment load using MCMC with optimal state number

Yang Zihan, Song Zhenghe※

(,,100083,)

The selection of state number depends highly on the subjective experience in the traditional Markov Chain Monte Carlo (MCMC). However, an inappropriate value of state number can lead to a great reduction in the accuracy of load simulation, even an increase in the running time during the simulation of agricultural equipment loads. This study aims to clarify the effect of state number on the simulation when the MCMC was applied to agricultural equipment load. Specifically, the mean error, standard deviation error, and deviation of rain flow matrix between the simulated and original load decreased rapidly to stabilize, as the state number increased. Moreover, the indicators were not generalizable, if there was no significance between them. An optimization of state number was also proposed using pseudo damage consistency. As such, the damage consistency between the simulated and original load gradually improved and smoothed out, as the state number increased, whereas, the rate of increase in the operation time continued to increase. The optimal state number was calculated to satisfy the damage consistency and minimum operation time, where a threshold value was set for the pseudo damage factor. Furthermore, the field tests were carried out for both tractor ploughing and soil preparation. The specific parameters were measured to validate, including the front axle vibration, front axle stress, and driveshaft torque load. The vibration loads were also utilized to apply for the tractor front drive axle during ploughing operations. It was found that the MCMC using optimal state number can be expected torealize the load simulation with pseudo damage differences within 1%. Furthermore, there were more significant differences between the load segments in the adjustment stage, where the optimal state numbers for each load segment were more dispersed than that in the operation stage. A cyclic simulation was also developed for the loads of key components, according to the operational characteristics of a tractor. Subsequently, the MCMC cycle simulations were also performed on the front axle vibration loads for ploughing. The results show that the simulated load retained the alternating switching between the operating and adjustment stages under tractor ploughing. The same procedure was used to simulate the stress load on the front axle under ploughing, where the torque was separately loaded on the driveshaft under soil preparation. The statistical characteristic indicators were selected, including the mean, standard deviation, and the maximum load cycle amplitude for each load segment. The deviation range of each statistical eigen value was also obtained, compared with the original. The eigen values simulation for each load segment was in a higher agreement with the original eigen values. The generality was further validated when applied to the load simulation of agricultural equipment with the objective of load spectrum preparation. Consequently, the MCMC using optimal state number was better matched to the target requirements of load spectrum preparation, compared with the conventional. The finding can also effectively reduce the computational cost for the higher accuracy during load simulation of agricultural machinery.

agricultural machinery; simulation; load; Markov chain; Monte Carlo method; optimal state number; pseudo damage

2021-08-20

2021-09-30

國(guó)家重點(diǎn)研發(fā)計(jì)劃資助項(xiàng)目（2017YFD0700301）

楊子涵，博士生，研究方向?yàn)檗r(nóng)機(jī)裝備載荷測(cè)試，載荷譜編制關(guān)鍵技術(shù)。Email：yangzihan@cau.edu.cn

宋正河，博士，教授，研究方向?yàn)檗r(nóng)機(jī)裝備試驗(yàn)驗(yàn)證方法與技術(shù)，智能化設(shè)計(jì)。Email：songzhenghe@cau.edu.cn

10.11975/j.issn.1002-6819.2021.20.002

S220

A

1002-6819(2021)-20-0015-08

楊子涵，宋正河. 利用優(yōu)選狀態(tài)數(shù)的MCMC 模擬農(nóng)機(jī)裝備負(fù)載[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2021，37(20)：15-22. doi：10.11975/j.issn.1002-6819.2021.20.002 http://www.tcsae.org

Yang Zihan, Song Zhenghe. Simulation of agricultural equipment load using MCMC with optimal state number[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(20): 15-22. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2021.20.002 http://www.tcsae.org