基于FDR閾值自動選取的拖拉機PTO轉(zhuǎn)矩載荷譜外推

王 禹 王 玲 呂東曉 溫昌凱 王書茂

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

0 引言

拖拉機作業(yè)環(huán)境復(fù)雜、載荷多變，這對其整機及零部件的疲勞壽命和可靠性造成較大影響。載荷譜作為研究作業(yè)載荷與疲勞壽命及可靠性關(guān)系的主要手段，能夠反映整機或零部件所承受的載荷時間歷程，是拖拉機進行疲勞及可靠性試驗的基礎(chǔ)。拖拉機動力輸出軸(Power take-off, PTO)是進行旋耕、耙地等作業(yè)時的主要動力輸出部件，其疲勞壽命及可靠性直接影響作業(yè)效率[1]。因此，編制適用于拖拉機PTO部件的動態(tài)載荷譜具有實際應(yīng)用價值。

由于受時間和測試成本的限制，一般多采用載荷譜外推方法來獲取全生命周期內(nèi)的載荷譜。常用載荷譜外推方法主要分為頻域外推和時域外推[2]。頻域外推出現(xiàn)時間較早，方法較為成熟[3]。如葛宜元等[4]、張英爽等[5]、李雯等[6]分別對整株秸稈還田機刀軸、工程車輛、汽車電機驅(qū)動系統(tǒng)的載荷譜編制方法進行了研究。頻域外推法需要將時域內(nèi)的載荷信號轉(zhuǎn)換到頻域，只能保留載荷循環(huán)的均幅值信息，會打破原有載荷的加載時序，對后續(xù)時域載荷譜應(yīng)用帶來不便。時域外推方法可以通過篩選對載荷信號擬合、直接外推，保留了載荷時間歷程，降低了轉(zhuǎn)換帶來的誤差。目前，常用時域外推方法主要基于極值理論，通過擬合超越閾值部分的數(shù)據(jù)來加以實現(xiàn)。JOHANNESSON[7]率先提出了基于極值理論的時域外推方法，并與雨流域外推的結(jié)果進行對比分析，驗證了其適用性與準確性。閾值的選取是時域外推準確性的關(guān)鍵，常用的閾值選取方法主要分為圖像法和數(shù)值法。圖像法多采用平均超出量函數(shù)圖[8]或Hill圖法[9]，但該方法具有較強的主觀性，閾值選取多憑經(jīng)驗確定，影響了擬合效果[10]。楊子涵等[11]基于超閾值模型對拖拉機傳動軸載荷進行了外推，采用圖像法確定候選閾值范圍，通過灰色關(guān)聯(lián)度分析方法對載荷數(shù)據(jù)中上下限的最優(yōu)閾值進行選取，并實現(xiàn)了外推。數(shù)值法則是基于統(tǒng)計學(xué)方法通過計算和分析選取合適的閾值，客觀性較好。如BEIRLANT等[12]用最小均方根誤差對應(yīng)的值作為閾值； THOMPSON等[13]引入一個新的參量，通過檢驗其分布特性實現(xiàn)了閾值自動選取，并進行了驗證；翟新婷等[14]從載荷數(shù)據(jù)擬合角度出發(fā)，將熵法、多準則決策技術(shù)結(jié)合，提出了基于多準則的閾值選取方法，并用于裝載機工作載荷譜的時域外推中。雖然相關(guān)學(xué)者對閾值自動選取方法進行了探討，但依舊存在計算精度較低、適用范圍有限及過程復(fù)雜的問題。

針對上述問題，本文基于錯誤發(fā)現(xiàn)率(False discovery rate, FDR)提出一種閾值自動選取方法，將多重假設(shè)檢驗應(yīng)用于閾值選取過程中，以實現(xiàn)閾值的自動快速選取。在此基礎(chǔ)上，以拖拉機PTO為研究對象進行轉(zhuǎn)矩載荷譜外推，并基于實際田間作業(yè)數(shù)據(jù)獲得拖拉機PTO全生命周期內(nèi)的轉(zhuǎn)矩載荷譜，通過統(tǒng)計學(xué)特征與雨流計數(shù)分析對閾值自動選取方法的準確性和時域外推方法的可行性進行驗證，以期為農(nóng)業(yè)機械作業(yè)載荷譜的外推研究提供參考。

1 基于FDR的閾值自動選取方法

本文基于超閾值模型實現(xiàn)拖拉機PTO轉(zhuǎn)矩載荷的時域外推，其原理是采用廣義帕累托分布(Generalized Pareto distribution，GPD)來擬合樣本中超越閾值數(shù)據(jù)的分布，根據(jù)擬合得到的GPD分布函數(shù)產(chǎn)生隨機序列對原始數(shù)據(jù)進行替換進而實現(xiàn)載荷外推，該方法可基于當前工況的極值信息外推可能存在的理論極值載荷，得到全生命周期內(nèi)的可能的轉(zhuǎn)矩及其頻次，實現(xiàn)極大值和極小值的雙向外推。

廣義帕累托分布是實現(xiàn)超閾值數(shù)據(jù)外推的基礎(chǔ)，其累積分布函數(shù)和概率密度函數(shù)表達式分別為

(1)

(2)

式中u——閾值x——樣本值

ξ——形狀參數(shù)

σ——尺度參數(shù)

GPD形狀參數(shù)和尺度參數(shù)的計算可以采用極大似然估計實現(xiàn)，但為了建立準確的超閾值模型，需要選取合適的閾值。閾值u太大，超出量個數(shù)少，導(dǎo)致有效樣本載荷數(shù)據(jù)丟失，使得GPD參數(shù)估計的方差變大；閾值u太小，雖然大量數(shù)據(jù)被保留下來，但是超出量數(shù)據(jù)與GPD分布相差較大，導(dǎo)致GPD參數(shù)估計偏差較大。為了提高載荷樣本擬合的準確性，本文將多重檢驗方法引入閾值選取過程，提出了基于FDR的閾值自動選取方法，實現(xiàn)了閾值自動快速選取。

通常，為了確定某閾值下的GPD分布擬合結(jié)果是否能夠準確反映樣本的分布特征，會采用假設(shè)檢驗的方法來進行判斷。P值是用來判定假設(shè)檢驗結(jié)果的參數(shù)，能夠反映在某一置信水平下零假設(shè)被接受的概率，可以理解為對某次假設(shè)檢驗可信程度的評估而不是保證結(jié)果為真，因此這意味著單次檢驗會存在一個缺陷：當只對單個閾值對應(yīng)的擬合結(jié)果進行假設(shè)檢驗時，如果假設(shè)檢驗的接受標準不具有很強的約束力，那么對于零假設(shè)的接受就可能會發(fā)生在較低的閾值，從而引入假陽性錯誤，導(dǎo)致閾值選取不合理[15]。FDR能夠表示在所有R次拒絕中錯誤發(fā)現(xiàn)的期望比例，因此BENJAMINI等[16]提出了FDR錯誤控制法(也稱為BH法)，基于多重檢驗方法，通過控制樣本整體零假設(shè)檢驗的FDR來決定P值的閾值，進而降低因檢驗次數(shù)增加而帶來的Ⅰ型錯誤率，但BH法會受制于樣本的檢驗順序?qū)е掠嬎懔鞒虖?fù)雜且約束力較弱[15]。G’SELL等[17]提出的ForwardStop規(guī)則采用順序多重檢驗的方法實現(xiàn)了FDR的控制，在保證檢驗準確性的基礎(chǔ)上提升了計算效率，因此，本文基于該規(guī)則并結(jié)合文獻[18]提出的擬合優(yōu)度檢驗方法實現(xiàn)閾值的自動選取，流程如下。

單次假設(shè)檢驗的P值計算采用了Anderson-Darling檢驗方法(AD檢驗)，該檢驗方法是二次經(jīng)驗分布函數(shù)統(tǒng)計的一種，可以用來檢驗數(shù)據(jù)與特定分布的擬合程度，引入了權(quán)重函數(shù)使得在計算過程中能夠更注重尾部分布的數(shù)據(jù)，對于檢測背離樣本的數(shù)據(jù)更加敏感，與其他檢驗方法(卡方檢驗、K-S檢驗)相比，更適用于極值樣本的檢驗[18]，計算式為

(3)

其中

w(x)=[F(x)(1-F(x))]-1

式中D——假設(shè)分布F(x)和經(jīng)驗分布函數(shù)Fn(x)之間的距離

n——樣本數(shù)w(x)——權(quán)重函數(shù)

(4)

其中

式中Z(i)——標準正態(tài)分布的累積分布函數(shù)

(5)

式中α——預(yù)期閾值，可取0.05、0.1、0.2，預(yù)期閾值越小，表示檢驗效果越嚴格，本文取0.05

依次對前i次的P值進行校正，并同時構(gòu)建校正后的P值序列；最后找到滿足預(yù)期閾值α最小的值，此時序列所對應(yīng)的閾值即為最優(yōu)閾值。

2 拖拉機PTO轉(zhuǎn)矩載荷超閾值模型建立

2.1 載荷譜采集試驗與數(shù)據(jù)預(yù)處理

2018年11月8—12日，在北京市郊區(qū)壤土地塊進行PTO載荷譜采集試驗。以TS404型拖拉機(山東五征集團有限公司)為樣機，試驗工況為旋耕，旋耕機型號為1GLN-0145，耕深為200 mm，耕作寬幅為1.8 m，平均作業(yè)速度為7.3 km/h，基于轉(zhuǎn)矩無線監(jiān)測系統(tǒng)對拖拉機PTO轉(zhuǎn)矩進行實時記錄，田間試驗現(xiàn)場如圖1所示。

轉(zhuǎn)矩無線監(jiān)測系統(tǒng)主要由無線雙法蘭扭矩傳感器、法蘭盤連接件、信號接收裝置、采集模塊和上位機軟件5部分組成。無線雙法蘭轉(zhuǎn)矩傳感器型號為CYB-807W，測試精度為10 N·m，該傳感器結(jié)構(gòu)緊湊，能夠與拖拉機PTO及萬向節(jié)傳動軸直接連接，并通過無線模塊輸出扭矩數(shù)據(jù)。轉(zhuǎn)矩信號通過無線接收器傳輸?shù)叫盘柌杉KNI9203，最終輸入到基于LabVIEW開發(fā)的上位機軟件平臺實時顯示與儲存，根據(jù)旋耕機刀具與土壤接觸頻率設(shè)置轉(zhuǎn)矩采樣頻率為20 Hz。轉(zhuǎn)矩無線監(jiān)測系統(tǒng)原理圖如圖2所示。

由于試驗環(huán)境復(fù)雜，PTO轉(zhuǎn)矩數(shù)據(jù)不可避免地會存在異常點，需要對數(shù)據(jù)進行預(yù)處理。選取作業(yè)工況良好的某次完整時域信號為外推樣本，利用Matlab對數(shù)據(jù)進行濾波去噪和去除趨勢項，并采用幅值門限法對數(shù)據(jù)進行預(yù)處理，剔除轉(zhuǎn)矩信號中的異常數(shù)據(jù)，提高載荷譜數(shù)據(jù)的可靠性。此外，實際采集到的載荷信號中還包含了許多對疲勞壽命影響很小的無效載荷循環(huán)，根據(jù)文獻[19]，一般選取最大循環(huán)幅值的10%作為基準對數(shù)據(jù)進行剔除，以簡化載荷時間歷程。圖3為處理后的拖拉機PTO轉(zhuǎn)矩載荷時間歷程，表1為載荷數(shù)據(jù)剔除前與剔除后的統(tǒng)計特性，可以看出，通過剔除無效載荷，數(shù)據(jù)量減少了50.14%，但其統(tǒng)計特性基本維持不變。

表1 剔除異常數(shù)據(jù)前后的PTO轉(zhuǎn)矩載荷數(shù)據(jù)統(tǒng)計特性Tab.1 Statistical property values of PTO torque load value

2.2 門限閾值確定與GPD擬合

采用自動閾值選取方法對處理后的PTO載荷數(shù)據(jù)進行門限閾值確定。對于載荷樣本中的極大值，選取候選閾值區(qū)間為[41 N·m，504 N·m]，以1 N·m為間隔遞增，基于上述流程進行自動計算，選取ForwardStop自動停止的序列代表的上限閾值為342 N·m，對應(yīng)超出量個數(shù)為640；同理，對載荷樣本中的下限閾值進行自動計算，候選閾值區(qū)間為[41 N·m，504 N·m]，以1 N·m為間隔遞減，選取ForwardStop自動停止的序列代表的下限閾值為100 N·m，超出量個數(shù)為159。超出上、下限閾值的樣本提取結(jié)果如圖4所示。

根據(jù)閾值選取結(jié)果，采用Matlab對超閾值樣本實現(xiàn)GPD參數(shù)的極大似然估計，得到上限閾值與下限閾值所對應(yīng)的GPD分布參數(shù)的估計值，如表2所示。

表2 GPD擬合結(jié)果 Tab.2 GPD fitting result

因此，可確定PTO載荷樣本對于上限閾值和下限閾值超出量數(shù)據(jù)的GPD分布函數(shù)分別為

(6)

(7)

2.3 不同閾值選取方法GPD擬合結(jié)果檢驗與對比

為進一步驗證和分析上述自動選取閾值方法的可行性以及擬合結(jié)果的精確性，同時采用平均超出量函數(shù)圖法對載荷數(shù)據(jù)進行閾值選取及GPD擬合，對兩者的擬合結(jié)果進行對比。

對于載荷樣本數(shù)據(jù)中的上限閾值，其樣本平均超出量如圖5所示。從圖中可以看到，由于超閾值樣本數(shù)量減少，樣本的平均超出量會隨著閾值的增大而劇烈波動，在閾值選取過程中需要盡量避免。根據(jù)曲線形狀，選取最接近波動區(qū)間且曲線分布近似為直線的起點作為最優(yōu)閾值，u=350 N·m。同理，對于載荷樣本中的下限閾值，確定其候選閾值為u=123 N·m。

針對每個閾值采用極大似然估計對超閾值數(shù)據(jù)進行參數(shù)估計，得到各自擬合后的GPD分布函數(shù)，從擬合檢驗、參數(shù)估計標準差、決定系數(shù)、擬合優(yōu)度4方面與圖像法選取的閾值進行量化分析與對比，結(jié)果如表3所示。

由表3可知，各閾值對應(yīng)的擬合結(jié)果在95%置

表3 不同閾值選取方法擬合結(jié)果對比Tab.3 Comparison of fitting results of different threshold selection methods

信水平下都通過了K-S檢驗和卡方檢驗，表明擬合結(jié)果準確，但自動閾值選取法所對應(yīng)的檢驗量高于圖像法，擬合結(jié)果較優(yōu)；通過參數(shù)擬合標準差來看，對于上限閾值，自動閾值選取法的參數(shù)擬合標準差更小，其擬合精度更高；而對于下限閾值，圖像法參數(shù)擬合標準差更小，但總體來看參數(shù)擬合標準差均很小，擬合效果可行；其次，通過決定系數(shù)可以看出，兩種方法的R2均大于0.995，表明擬合結(jié)果與超閾值樣本具有較高的重合度，效果良好；從擬合優(yōu)度來看，對于上限閾值，自動閾值選取法的擬合優(yōu)度比圖像法的擬合優(yōu)度小8.7%，而對于下限閾值，自動閾值選取法的擬合優(yōu)度比圖像法的擬合優(yōu)度小31.21%，表明自動閾值選取法的擬合結(jié)果在總體上更接近樣本數(shù)據(jù)。通過以上對比，綜合各項數(shù)據(jù)來看，基于FDR的自動閾值選取法優(yōu)于傳統(tǒng)的圖像法。

為了更加直觀反映自動閾值選取法得到的擬合效果，繪制樣本載荷與擬合得到的GPD分布函數(shù)的累積分布函數(shù)圖(Cumulative distribution function，CDF)和Q-Q圖進行初步驗證。如圖6所示，兩者CDF圖具有很高的重合度，相關(guān)系數(shù)大于0.99，同時Q-Q圖的分布也近似為一條直線，從側(cè)面表明了擬合結(jié)果的準確性。

3 全生命周期內(nèi)PTO轉(zhuǎn)矩載荷譜外推與驗證

3.1 載荷譜外推

根據(jù)式(6)和式(7)，利用Matlab生成隨機數(shù)據(jù)，并在原位置對樣本數(shù)據(jù)進行替換得到外推因子為1的時域外推載荷譜，原載荷譜與1倍外推載荷譜對比如圖7所示。圖8是從1倍外推載荷譜中截取的某一段載荷示意圖，可以看出，原載荷樣本中超過上限閾值和低于下限閾值的數(shù)據(jù)被GPD擬合生成的隨機數(shù)據(jù)替換，而處于上限閾值與下限閾值區(qū)間內(nèi)的數(shù)據(jù)則被保留下來，且新生成的極值數(shù)據(jù)有的大于原載荷數(shù)據(jù)，有的小于或等于原載荷數(shù)據(jù)，表明采用時域外推方法可以有效保留原有的載荷時序，而且還可以得到載荷時間歷程中未出現(xiàn)過的極值載荷，實現(xiàn)載荷數(shù)據(jù)的外推。

為了得到所有可能出現(xiàn)的載荷數(shù)據(jù)，需要將載荷累積頻次外推至106次。重復(fù)上述步驟，產(chǎn)生隨機數(shù)對原載荷數(shù)據(jù)進行外推，分別得到外推因子為1、50、131的時域外推數(shù)據(jù)，各外推載荷數(shù)據(jù)的累積頻次如圖9所示，當外推因子為131時，載荷的累積頻次外推至106次。

3.2 外推載荷譜驗證

為了驗證時域外推載荷結(jié)果的合理性，對外推載荷譜進行驗證。

3.2.1統(tǒng)計特性分析

對外推得到的載荷數(shù)據(jù)進行統(tǒng)計特性分析，如表4所示。從表中可以看出，隨著外推數(shù)據(jù)量的增加，載荷極值范圍從[40.01 N·m, 504.20 N·m]擴大到了[36.12 N·m, 539.19 N·m]，載荷覆蓋范圍擴大了8.36%；同時可以看出，外推載荷數(shù)據(jù)的均值、標準差以及方差基本保持了原載荷樣本的統(tǒng)計特性，最大限度地保留了原有的載荷特征，可以認為外推載荷譜是有效的。

表4 外推后的PTO轉(zhuǎn)矩載荷數(shù)據(jù)統(tǒng)計特性Tab.4 Statistical property values of PTO torque load value after extrapolation

3.2.2雨流計數(shù)分析

應(yīng)用雨流計數(shù)法對原載荷數(shù)據(jù)、1倍外推載荷數(shù)據(jù)以及最終載荷數(shù)據(jù)進行均幅值雙參數(shù)統(tǒng)計。原載荷與1倍外推載荷數(shù)據(jù)的均值、幅值頻次對比如圖10所示，對二者的分布進行相關(guān)性分析，得到幅值分布的相關(guān)系數(shù)為0.995，均值分布的相關(guān)系數(shù)為0.989，可見1倍外推載荷數(shù)據(jù)與原載荷數(shù)據(jù)的分布具有極大的相似性，較好地模擬了原載荷的數(shù)據(jù)分布特點及規(guī)律。

由于最終外推載荷數(shù)據(jù)量與原載荷數(shù)據(jù)量不同，因此基于雨流矩陣分布情況對外推結(jié)果進行分析。如圖11所示，相對于原載荷數(shù)據(jù)，可以看到對原載荷循環(huán)的均值和幅值以及頻次實現(xiàn)了同步外推。如圖11中紅色圈所示，經(jīng)過外推后的載荷數(shù)據(jù)在極大值與極小值兩端數(shù)據(jù)的分布出現(xiàn)了2個山峰，意味著將載荷數(shù)據(jù)的均幅值及對應(yīng)的頻次實現(xiàn)了擴充；而從雨流矩陣的整體分布來看，外推后的載荷譜與原載荷樣本的雨流矩陣具有一致的分布規(guī)律，證明該時域外推方法是可行合理的。

4 結(jié)論

(1)基于旋耕作業(yè)工況采集拖拉機PTO實際作業(yè)載荷，并利用開發(fā)的PTO轉(zhuǎn)矩無線監(jiān)測系統(tǒng)得到了拖拉機實際作業(yè)載荷數(shù)據(jù)，通過對數(shù)據(jù)進行預(yù)處理和剔除無效載荷，提高了載荷數(shù)據(jù)的準確性和利用率。

(2)針對時域外推過程閾值選取主觀性強的問題，基于FDR將AD檢驗與多重檢驗相結(jié)合，實現(xiàn)了閾值的自動選取，確定了拖拉機PTO扭矩時域外推的上限閾值為342 N·m，下限閾值為100 N·m。從擬合檢驗、標準差、決定系數(shù)、擬合優(yōu)度4方面與圖像法選取的閾值進行了對比與檢驗，結(jié)果表明，兩種閾值的擬合結(jié)果在95%置信水平下都通過了K-S檢驗和卡方檢驗，且參數(shù)擬合的標準差都比較小，總體來說自動閾值選取法的擬合結(jié)果相對較優(yōu)；兩種方法的擬合結(jié)果的決定系數(shù)都比較高(R2>0.995)，表明均與樣本有較高的重合度；從擬合優(yōu)度來看，對于上限閾值，自動閾值選取法的擬合優(yōu)度比圖像法的擬合優(yōu)度降低8.7%，而對于下限閾值，自動閾值選取法的擬合優(yōu)度比圖像法的擬合優(yōu)度降低31.21%，表明自動閾值選取法的擬合結(jié)果在總體上更接近樣本數(shù)據(jù)。

(3)基于超閾值模型實現(xiàn)了拖拉機PTO轉(zhuǎn)矩載荷譜的時域外推。當外推因子為131、轉(zhuǎn)矩累積頻次達到106次時，可認為其能夠反映全生命周期內(nèi)的載荷歷程；從載荷數(shù)據(jù)的統(tǒng)計特性和雨流計數(shù)分析兩方面對外推載荷譜的可行性進行了分析與驗證。外推后的載荷譜在統(tǒng)計特性方面與原有載荷譜相比變化較小，在實現(xiàn)極值外推的過程中最大限度地保留了原有的載荷特征，可認為是有效的；在雨流計數(shù)分析中，原載荷與1倍外推載荷頻次分布的均值、幅值的相關(guān)系數(shù)分別為0.989和0.995，說明了時域外推方法的準確性。外推載荷譜與原載荷樣本雨流矩陣具有相似的分布規(guī)律，可以證明該時域外推方法合理可行。