基于DTW距離的拖拉機(jī)傳動(dòng)軸載荷樣本長(zhǎng)度計(jì)算方法

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

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

0 引 言

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

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

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

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

1 基于DTW距離的載荷樣本長(zhǎng)度計(jì)算

1.1 DTW距離基本原理

DTW距離是將時(shí)間序列拉伸或收縮來進(jìn)行對(duì)應(yīng)相似點(diǎn)距離的累加[22-23]，是評(píng)價(jià)時(shí)間序列特征相似性的一種方法，基本原理如下。

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

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

1.2 計(jì)算流程

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

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

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

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

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

2 拖拉機(jī)傳動(dòng)軸載荷數(shù)據(jù)測(cè)試及分析

2.1 扭矩信號(hào)獲取

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

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

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

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

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

圖4 拖拉機(jī)傳動(dòng)軸犁耕作業(yè)扭矩信號(hào)Fig.4 Torque signal of tractor drive shaft during ploughing operation

2.2 扭矩信號(hào)分析

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

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

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

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

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

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

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

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

3 實(shí)例分析

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

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

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

3.1 近似均值精度估計(jì)法確定載荷樣本長(zhǎng)度

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

式中S(x)為載荷樣本標(biāo)準(zhǔn)差，MPa；為載荷樣本均值，MPa；N為相互獨(dú)立的子樣本個(gè)數(shù)；x為不同子樣本個(gè)數(shù)的載荷均值，MPa。

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

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

3.2 均值曲線擬合法確定載荷樣本長(zhǎng)度

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

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

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

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

3.3 基于DTW距離確定傳動(dòng)軸載荷樣本長(zhǎng)度

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

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

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

3.4 載荷樣本長(zhǎng)度驗(yàn)證

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

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

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

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

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

表2 載荷樣本長(zhǎng)度計(jì)算方法對(duì)比Table 2 Comparison of load sample size calculation methods

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

4 結(jié) 論

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

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

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