尹宜勇，付寧善，廖 頻，宋正河

（中國農業(yè)大學現(xiàn)代農業(yè)裝備優(yōu)化設計北京市重點實驗室，北京 100083）

0 引 言

載荷譜編制是疲勞壽命分析和疲勞可靠性試驗的關鍵環(huán)節(jié)[1-3]，其中載荷樣本長度的確定是編制載荷譜中的關鍵步驟[4-5]。拖拉機田間工作時，傳動軸向轉向驅動橋傳遞動力，其是影響拖拉機整機性能的重要因素[6]。因此，以拖拉機傳動軸為研究對象，確定載荷樣本長度，對編制拖拉機傳動軸田間作業(yè)工況下的載荷譜具有重要意義。

工程領域確定載荷樣本長度的方法主要有近似均值精度估計法、曲線擬合法[7-8]、疲勞壽命法[9-10]以及基于貝葉斯算法的載荷樣本長度計算方法[11]。由于農業(yè)機械作業(yè)環(huán)境復雜開放[12]以及地塊大小的影響，農業(yè)機械所受載荷波動較大且工作階段占比會發(fā)生變化，使載荷的統(tǒng)計特征發(fā)生較大變化，而近似均值精度估計法和曲線擬合法主要依靠數(shù)據(jù)統(tǒng)計特征進行載荷樣本長度的計算，因此會有較大誤差；疲勞壽命法的計算結果會因疲勞理論與載荷本身特點的不同而產(chǎn)生較大差異[13]；載荷均值服從混合分布[14]，使貝葉斯算法計算精度較低。

載荷譜本質上是反應整機結構或關鍵零部件受載情況的載荷時間歷程[15-16]，而確定載荷樣本長度是選取能夠代表總體載荷特征的最短載荷時間序列，所以確定載荷樣本長度可認為是對載荷時間序列相似性程度的度量。動態(tài)時間扭曲（DTW, Dynamic Time Warping）距離主要應用于語音識別和在線簽名驗證領域[17-19]，是對于時間序列數(shù)據(jù)特征相似性度量的一種方法[20]，其能夠進行不等長時間序列相似性判別、運用特征匹配進行相似性判斷[21]，能夠確定載荷樣本與載荷總體相似性程度且受載荷統(tǒng)計特征的影響較小。

鑒于DTW距離方法的特點，本文針對工程領域載荷樣本長度計算方法在農業(yè)機械領域存在的問題，提出一種基于DTW距離計算載荷樣本長度的方法。以拖拉機傳動軸載荷為研究對象，獲得拖拉機犁耕時傳動軸的田間動態(tài)載荷數(shù)據(jù)，并根據(jù)犁耕作業(yè)特點進行作業(yè)工況劃分，運用DTW距離計算載荷樣本長度，并運用參數(shù)外推中的均值服從混合正態(tài)分布的擬合參數(shù)相對誤差進行檢驗，驗證基于DTW距離載荷樣本長度計算方法的適用性。

1 基于DTW距離的載荷樣本長度計算

1.1 DTW距離基本原理

DTW距離是將時間序列拉伸或收縮來進行對應相似點距離的累加[22-23]，是評價時間序列特征相似性的一種方法，基本原理如下。

選取2個時間序列S1(i)、S2(j)，長度分別為r、s，由這2個時間序列建一個r×s的矩陣，則這2個數(shù)據(jù)序列上任意兩點之間的DTW距離為：

式中d(i,j)為序列點S1(i)和S2(j)間的距離（一般為歐氏距離 ） ； γ (i- 1 , j - 1 )為 從 元 素 (S1( 1),S2(1))到 元 素(S1(i- 1 ),S2(j- 1 ))間的最小累計距離；γ(i- 1 ,j)為從元素(S1( 1),S2(1))到元素 (S1(i- 1 ),S2(j))間的最小累計距離；γ(i,j- 1 )為從元素 (S1( 1),S2(1))到元素 (S1(i),S2(j- 1 ))間的最小累計距離；γ(i,j)為從元素 (S1( 1),S2(1))到元素(S1(i),S2(j))間的最小累計距離。

1.2 計算流程

圖1 基于DTW距離確定載荷樣本長度流程圖Fig.1 Flow chart for determining the load sample size with Dynamic Time Warping (DTW) distance

DTW距離是運用動態(tài)時間規(guī)劃的方法，將時間序列特征進行匹配，從而進行相似性判別，不參照數(shù)據(jù)的統(tǒng)計特征，因此，針對農業(yè)機械的載荷樣本長度計算較為適用。應用基于DTW方法進行載荷樣本長度的計算流程如圖1所示。隨著子樣本個數(shù)的不斷疊加，載荷時間序列中樣本點的個數(shù)會增加，會造成兩載荷時間序列之間的DTW距離增大，為保證不同子樣本個數(shù)載荷時間序列的DTW距離處于同一評價標準下，按照式（2）對不同樣本點數(shù)量的載荷時間序列進行歸一化處理。

式中S(i)為原時間序列；Smax為原時間序列當中的最大值；Smin為原時間序列當中的最小值；Ns為原時間序列中樣本點數(shù)量；S′(i)為歸一化后的時間序列。

將x個子樣本經(jīng)歸一化處理之后得到的樣本載荷序列 （x1,… ,xi, … ,xm） 與 歸 一 化 后 的 總 體 載 荷 序 列（ y1, … , yj, … , yn）數(shù)據(jù)點構建一個n×m距陣，之后計算樣本載荷數(shù)據(jù)點xi與總體載荷數(shù)據(jù)點yj之間的歐氏距離，并放入矩陣中點（i,j）處。然后根據(jù)式（1）計算樣本載荷序列與總體載荷序列之間的DTW距離。計算時應該注意 DTW 距離的計算應該從矩陣單元的成對角的起始點處開始和終結點處結束，并保證經(jīng)過的路徑連續(xù)、單調，使樣本載荷序列和總體載荷序列的特征相對應，獲得兩者之間的最短距離。

得到距離之后，選取給定的誤差εr，將子樣本不斷累加并與載荷總體y進行DTW距離d的計算，若子樣本累加到使d小于εr，則輸出子樣本數(shù)量，從而確定出合適的載荷樣本長度。

2 拖拉機傳動軸載荷數(shù)據(jù)測試及分析

2.1 扭矩信號獲取

由于田間道路的路面不平度較高，拖拉機整機振動明顯，而且傳動軸需要高速旋轉，因此，扭矩采集系統(tǒng)需保證方便安裝并具有一定的抗振性能。本試驗采用北京必創(chuàng)公司生產(chǎn)的TQ201型無線扭矩節(jié)點進行扭矩信號測取，其具有一定抗振性能，且體積較小、安裝方便。將無線扭矩節(jié)點、BF350-3BA型半橋應變片、配套無線接收器、供電單元以及筆記本共同構成無線扭矩采集系統(tǒng)采集傳動軸扭矩，其中無線扭矩節(jié)點布置如圖2所示。

圖2 無線扭矩節(jié)點布置Fig.2 Arrangement of wireless torque node

試驗時間為2018年10月，試驗地點為河南省洛陽市孟津縣金村，選取面積約13.3 hm2，長度約為170 m，玉米收獲完成后的農田為試驗場地進行犁耕作業(yè)，土壤表面玉米秸稈留茬高度約為10 cm，土壤種類為沙土。按照犁耕標準[24]對拖拉機犁耕工況下的作業(yè)環(huán)境和作業(yè)質量進行檢測，測得環(huán)境溫度為25.6 ℃，環(huán)境濕度為21%，風速為4.6 m/s，犁耕作業(yè)幅寬為2.1 m，耕深為320 mm，犁耕后的碎土率為99.2%，經(jīng)檢驗作業(yè)質量符合國家標準。

犁耕過程中拖拉機擋位為中三擋，速度為 0～10 km/h，拖拉機田間犁耕作業(yè)如圖3所示。選取采樣頻率為100 Hz，在進行扭矩信號采集時，為了提高采集系統(tǒng)精度，采集過程中運用巴特沃斯濾波的方式對數(shù)據(jù)進行抗混疊濾波，降低環(huán)境帶來的干擾，故不需要再對載荷時間序列進行預處理。測取犁耕作業(yè)過程中拖拉機傳動軸扭矩的變化如圖4所示。

圖3 拖拉機田間犁耕作業(yè)Fig.3 Tractor ploughing in the field

圖4 拖拉機傳動軸犁耕作業(yè)扭矩信號Fig.4 Torque signal of tractor drive shaft during ploughing operation

2.2 扭矩信號分析

無線扭矩采集系統(tǒng)測取的載荷信號為扭矩信號，為方便后續(xù)分析處理，將扭矩載荷轉化為切應力載荷τ(MPa)。

式中 WP=/16為抗扭截面系數(shù)；da為傳動軸直徑，取30 mm；M為傳動軸所受扭矩，N? m。

為使試驗田犁耕完整，拖拉機在地塊邊界需進行姿態(tài)調整，由于此時變速箱向轉向驅動橋傳遞動力較小，且前進方向、擋位發(fā)生多次變化，所以此時傳動軸所受的切應力較小且有較大波動；在犁耕作業(yè)過程中，變速箱向轉向驅動橋傳遞動力較大，但拖拉機狀態(tài)相對平穩(wěn)，所以傳動軸承受的切應力較大且較平穩(wěn)。由于拖拉機在犁耕作業(yè)時和土地邊界調整時傳動軸扭矩載荷特征差異明顯，為使犁耕作業(yè)滿足各態(tài)歷經(jīng)性，使測量的載荷能夠代表田間作業(yè)情況下載荷總體特征[25]，將犁耕工況與調整工況兩種工況劃分為一個作業(yè)循環(huán)。拖拉機進行一個作業(yè)循環(huán)傳動軸受到的切應力如圖5所示。

為分析傳動軸扭矩頻率變化范圍，應用Matlab軟件進行頻譜分析，結果如圖 6所示。分析可知，拖拉機傳動軸受到的載荷主要是低頻信號，而由于農業(yè)機械作業(yè)的地面環(huán)境較為惡劣，傳動軸受到地面和變速箱的激勵，使得載荷的幅值在0.01和2 Hz附近有2個峰值。

圖5 作業(yè)工況劃分Fig.5 Division of operation conditions

圖6 載荷頻域分析Fig.6 Analysis of load frequency domain

在作業(yè)循環(huán)中調整工況占的比例較小，且會因為駕駛員駕駛習慣以及土地邊界的變化造成調整階段占比發(fā)生變化。選取前4個作業(yè)循環(huán)進行分析，結果如表1所示。調整工況相比犁耕工況載荷波動較大，這會使得每個作業(yè)循環(huán)的統(tǒng)計特性變化較大，造成傳統(tǒng)近似均值精度估計法和曲線擬合法對于拖拉機傳動軸載荷樣本長度的計算適用性較差。

表1 前4個作業(yè)循環(huán)統(tǒng)計分析結果Table 1 Statistical analysis results of the 4 cycles of work

3 實例分析

將一個作業(yè)循環(huán)作為一個子樣本，選取同一地段連續(xù)犁耕作業(yè) 20個子樣本長度的載荷時間序列為分析對象，時間長度為3 225 s。

由于小載荷循環(huán)對于零部件疲勞損傷的貢獻較小，一般可將低于最大載荷循環(huán) 10%的小載荷循環(huán)濾除[26]，使用nCode軟件將小載荷循環(huán)濾除。

只有載荷時間序列滿足平穩(wěn)性檢驗和各態(tài)歷經(jīng)性檢驗，才能用載荷樣本代替總體[25]。經(jīng)驗證，所測量犁耕狀態(tài)下傳動軸扭矩通過ADF（Augmented Dickey-Fuller）單位根檢驗，符合平穩(wěn)性，并且每個作業(yè)循環(huán)包括犁耕工況和調整工況，包含作業(yè)全過程，符合各態(tài)歷經(jīng)性要求，因此，能夠代表載荷總體進行載荷樣本長度的計算。

3.1 近似均值精度估計法確定載荷樣本長度

根據(jù)統(tǒng)計誤差的定義，再結合樣本長度與統(tǒng)計誤差εr的函數(shù)關系，選取載荷均值為樣本，得到置信度為95.4%下，估算載荷樣本長度的公式為[25]：

式中S(x)為載荷樣本標準差，MPa；為載荷樣本均值，MPa；N為相互獨立的子樣本個數(shù)；x為不同子樣本個數(shù)的載荷均值，MPa。

逐漸增加子樣本個數(shù)，得到不同子樣本個數(shù)的載荷均值如圖7所示。一般選取誤差εr為0.1[27],計算得到S(x)為6.511 MPa，為48.57 MPa，由式（4）得N=7.19，取N為8，即由近似均值精度估計法得到的最少子樣本數(shù)量為8。

圖7 載荷均值曲線擬合Fig.7 Fitting Curve of mean load

3.2 均值曲線擬合法確定載荷樣本長度

由圖 7可知，不同子樣本個數(shù)的載荷均值大致符合函數(shù)x=aNb+c，運用非線性函數(shù)x=f(N)擬合，在 95%置信水平下得到a=-27.66，b=-2.231，c=50.58，擬合優(yōu)度R2為0.9034，擬合效果良好。當子樣本個數(shù)N趨近于無窮時，載荷均值極限值x∞=x= 5 0.58 MPa 。根據(jù)中心極限定理，采樣容量N＞4，就可認為子樣均值抽樣分布接近于正態(tài)分布，取置信水平為 95%，得到均值估計的置信區(qū)間為

計算得到載荷均值的置信區(qū)間為[45.527，51.621]，x∞位于區(qū)間內，證明了極限x∞作為載荷均值總體參數(shù)的準確性。

獲得載荷均值總體參數(shù)后，同樣選取誤差εr=0.1，由式（6）可得載荷子樣本個數(shù)

經(jīng)計算得N為2.061，取N為3。即由載荷均值曲線擬合法得到的最少子樣本數(shù)量為3。

3.3 基于DTW距離確定傳動軸載荷樣本長度

將20個子樣本的載荷時間序列經(jīng)過歸一化處理作為總體。計算經(jīng)歸一化處理后不同子樣本個數(shù)的載荷序列與總體的DTW距離，得到結果如圖8所示。由于農田作業(yè)環(huán)境復雜開放，使拖拉機在進行犁耕作業(yè)時每個作業(yè)循環(huán)傳動軸載荷會有所差異，會出現(xiàn)如5個子樣本、6個子樣本與總體之間的 DTW 距離隨子樣本個數(shù)增多變大的情況。但總體來看，當子樣本個數(shù)不斷增加時，載荷時間序列與總體之間的DTW距離越來越小，表明樣本載荷與總體的相似性程度越來越高。

圖8 不同子樣本個數(shù)載荷序列與總體載荷序列的DTW距離Fig.8 DTW distance between load sequences of different number of sub-samples and total load sequence

為保證載荷樣本長度保持合適精度，并與傳統(tǒng)方法精度相同，選取DTW距離為0.1，由圖8可知，確定載荷樣本長度為16個子樣本。

3.4 載荷樣本長度驗證

目前針對載荷樣本長度的驗證方法主要是對載荷的統(tǒng)計特征進行驗證，沒有考慮載荷均幅值及其頻次，會對載荷外推產(chǎn)生影響，進而影響載荷譜編制精度。

在載荷外推中，通常采用參數(shù)外推進行載荷均幅值的頻次外推[28-29]，其實質是獲得載荷均幅值聯(lián)合分布函數(shù)，其中對載荷均值及其頻次采用混合正態(tài)分布進行擬合是關鍵步驟。因此，為保證載荷譜編制精度，本文以載荷均值及其頻次服從混合正態(tài)分布函數(shù)為依據(jù)[14]，進行載荷樣本長度的驗證。

將子樣本個數(shù)不斷疊加得到的載荷時間序列進行雨流計數(shù)法計數(shù)，然后對每個載荷時間序列的均值劃分為64級[30-31]。為簡化求解過程，降低參數(shù)復雜度，本文采用雙正態(tài)分布函數(shù)進行載荷均值及其頻次擬合，擬合函數(shù)如式（7）所示。

式中λ1、λ2為擬合系數(shù)；μ1、μ2為雙正態(tài)分布數(shù)學期望；σ1、σ2為雙正態(tài)分布標準差。

為保證擬合參數(shù)真實有效，在置信水平為 95%且保證擬合優(yōu)度大于0.9的前提下，采用Trust-Region方法將不同子樣本個數(shù)的載荷時間序列進行雙正態(tài)分布擬合，選取 20個子樣本載荷時間序列擬合之后的參數(shù)為參考值，計算由不同載荷樣本長度計算方法得到的載荷樣本長度對應的均值及其頻次雙正態(tài)分布參數(shù)的相對誤差，以此對不同載荷樣本長度計算方法進行對比（表2）。

表2 載荷樣本長度計算方法對比Table 2 Comparison of load sample size calculation methods

由表 2可知，由近似均值精度估計法和均值曲線擬合法得到的載荷樣本長度使得參數(shù)擬合外推中的均值及其頻次雙正態(tài)分布擬合參數(shù)相對誤差最大值分別為126.06%和80.62%，而基于DTW距離計算方法的擬合參數(shù)相對誤差最大值為5.43%，并能夠保證擬合參數(shù)相對誤差均小于10%，驗證了基于DTW距離的載荷樣本長度計算方法的適用性。

4 結 論

1）按照拖拉機傳動軸犁耕作業(yè)下載荷特性將拖拉機犁耕時作業(yè)工況和調整工況作為一個完整的犁耕作業(yè)循環(huán)，使載荷滿足各態(tài)歷經(jīng)性要求。

2）基于DTW距離進行載荷樣本長度的計算。為排除樣本點個數(shù)對于載荷距離的影響，將載荷進行歸一化，將載荷樣本不斷疊加的載荷時間序列與總體進行 DTW距離的計算，選取DTW距離為0.1，確定載荷樣本長度為16個子樣本。

3）基于參數(shù)外推當中均值及其頻次服從混合正態(tài)分布進行載荷樣本長度驗證。由近似均值精度估計法和均值曲線擬合法確定的載荷樣本長度的擬合參數(shù)相對誤差分別為126.06%和80.62%，而基于DTW距離方法確定的載荷樣本長度使擬合參數(shù)的相對誤差為5.43%，從而驗證了基于DTW距離的載荷樣本長度計算方法的適用性。