蘇力德 張 永 王 健 尹 玉 宗哲英 鞏彩麗

(1.內蒙古農業(yè)大學機電工程學院， 呼和浩特 010018； 2.內蒙古大學電子信息工程學院， 呼和浩特 010021)

0 引言

內蒙古自治區(qū)是我國重要的奶牛養(yǎng)殖優(yōu)勢區(qū)域，其奶牛存欄數、牛奶產量，以及乳制品企業(yè)的規(guī)模、效益和影響力均位于全國前列。據我國奶業(yè)年鑒統計，截止到2017年，奶牛存欄數已達到216萬頭，牛奶產量達到730萬噸，乳制品產量達到337萬噸。奶牛疾病阻礙了奶業(yè)的健康發(fā)展，據統計，我國大中型牧場臨床跛行率和嚴重跛行率分別為(31±12)%和(10±6)%[1]。跛行不僅導致奶牛身體狀況下降、牛奶產量降低，還將使奶牛過早地被淘汰。因此，準確分割奶牛步態(tài)、識別奶牛跛行對減少養(yǎng)殖場的經濟損失、提高福利化養(yǎng)殖水平具有重要意義。

步態(tài)由周期性的步幅組成，如何從步態(tài)序列中分割單個步幅已成為運動分析領域的研究熱點。國內外學者主要采用峰值檢測法[2-3]、視頻分割法[4]、模板匹配方法分割步態(tài)。峰值檢測法必須滿足一定的閾值條件才可以進行分割；視頻分割方法工作量大、自動化程度不高，限制了在實際中的應用；模板匹配方法通過計算兩個時間序列之間相似程度對步態(tài)進行分割。國內外基于模板匹配方法的步態(tài)分割主要集中在人體步態(tài)研究，對奶牛步態(tài)分割的研究尚未見報道。文獻[5-6]利用基于模板匹配的互相關方法對人體步態(tài)進行分割，但步幅模板長度和形式固定，不適用于不同步幅長度序列的分割。

動態(tài)時間規(guī)整算法(Dynamic time warping,DTW)是一種非線性匹配算法，該算法來源于動態(tài)規(guī)整原理，解決時間長度不一致時的模板匹配問題。文獻[7]利用動態(tài)時間規(guī)整算法對人體步態(tài)進行分割，通過研究步態(tài)周期數據間的差異性，證明該方法在臨床疾病檢測上優(yōu)于其他方法。文獻[8]利用傳感器技術采集步態(tài)信息，以動態(tài)時間規(guī)整算法分割步態(tài)，通過k近鄰分類器對偏癱患者進行識別，得出曲線下面積為0.94。文獻[9]利用動態(tài)時間規(guī)整算法對帕金森患者步態(tài)進行分割，獲得100%的分割準確率。文獻[10]利用動態(tài)時間規(guī)整算法對25名老年人的步態(tài)進行分割，分割準確率達到97%。

常規(guī)的動態(tài)時間規(guī)整算法樣本序列是一維數據，而奶牛步態(tài)加速度信號的時間和幅度是周期性變化的。本文利用三維加速度傳感器采集奶牛趾蹄加速度信號，將時間和幅度組成的二維數據作為步態(tài)序列，提出一種改進的動態(tài)時間規(guī)整算法，對奶牛步態(tài)進行分割，并對每個步幅提取步態(tài)特征，利用Wilcoxon秩和檢驗分析健康奶牛和跛行奶牛特征值之間的關系，選取與跛行相關的6個特征，利用邏輯回歸方法建立跛行識別模型，以準確識別奶牛的跛行。

1 材料和方法

1.1 測量設備

實驗場地為內蒙古旗幟牧業(yè)有限公司，選擇牧坤電子科技有限公司生產的M63X系列三維加速度傳感器，測量范圍±8g，采樣頻率50 Hz。將傳感器封裝在密閉防水盒，固定在制作好的綁帶上，形成可穿戴式設備。數據通過ZigBee協議無線傳輸到基站，并實時上傳到計算機。奶牛進入長、寬分別為23、0.8 m的測試通道，射頻識別閱讀器識別身份信息，使用PCO.dimax系列高速攝像機拍攝奶牛行走視頻，如圖1所示。由于奶牛85%的趾蹄疾病發(fā)生在后趾外側爪上[11]，因此將加速度傳感器綁在后肢外側，為避免奶牛應激反應產生的測量誤差，每頭牛綁上加速度傳感器行走10 min后開始采集，以此適應測量設備。

圖1 測試環(huán)境及設備Fig.1 Test environment and equipment1.加速度測量裝置 2.捆綁式加速度傳感器 3.無線接收基站

1.2 數據采集與預處理

數據采集使用實驗室虛擬儀器工程平臺(Laboratory virtual instrumentation engineering workbench，LabVIEW)實時顯示在上位機，LabVIEW與Access數據庫連接實現對數據的管理，數據處理及算法實現通過Matlab軟件編程，并與LabVIEW的DLL(Dynamic-link library)及ActiveX等技術實現混合編程。選取30頭荷斯坦奶牛，其中健康奶牛21頭、跛行奶牛9頭，在測量通道內連續(xù)采集6 d，每天上午和下午擠完奶各采集一次，每次重復采集一次。為確保數據的準確性高速攝像機與加速度傳感器系統實現同步采集，其中高速攝像機以940 f/s速率拍攝奶牛行走視頻，并以該視頻數據作為黃金標準。3名獸醫(yī)通過觀察每頭奶牛行走視頻，以5分制[12]評分規(guī)則對其進行評分，取平均值作為跛行評分值。在三維加速度傳感器中，定義X軸為垂直方向、Y軸為前進方向、Z軸為水平方向。

由于噪聲的存在，加速度數據不能直接應用于步態(tài)分割中，需對其進行預處理。小波變換去噪方法能有效區(qū)分信號突變部分和噪聲[13]，本文通過小波變換方法對加速度數據進行預處理。利用Matlab軟件在Daubechies小波基框架內，選擇DB(Daubechies)6小波進行3層分解，通過信號重構去除干擾信號。

2 改進動態(tài)時間規(guī)整算法步態(tài)分割實現

將連續(xù)的步態(tài)序列分割為單個步幅稱為步態(tài)分割。對于周期性步態(tài)的分割相對簡單實現，然而在現實中步幅并不完全是周期性的，不同跛行評分、不同體質量和泌乳的奶牛步幅行走模式均有所不同。這些異構的步態(tài)序列是步態(tài)分割的主要難點之一，而動態(tài)時間規(guī)整算法主要解決序列長度不一致時模板匹配的問題，適用于不同行走模式下的步態(tài)分割。奶牛步幅定義為同一肢體連續(xù)2次接觸地面之間的間隔，對應的時間差為一個步態(tài)周期。步態(tài)周期分為支撐階段和擺動階段，其中，在支撐階段，趾蹄與地面接觸開始，并在推開后終止。在擺動階段，趾蹄需在地面以上抬起。一個步態(tài)周期中，擺動時間占31.6%，支撐階段占68.4%[14]，平均步態(tài)周期為(1.22±0.05) s[15]。通過前期實驗分析可知，利用垂直方向加速度信號分割步態(tài)容易存在誤差，而水平加速度很大程度上受到乳房的影響，整個步行過程中變化并不明顯。前進方向加速度信號極大值點突出容易區(qū)分，本文采用該方向加速度信號對步態(tài)進行分割。

利用動態(tài)時間規(guī)整算法分割步態(tài)是通過求模板樣本序列映射到測試樣本序列的成本距離，利用遞推公式計算滿足最小累積成本的規(guī)整路徑，尋找與模板樣本相似的子序列，最終從連續(xù)的步態(tài)序列中分割單個步幅。本文算法的步態(tài)分割整個流程如圖2所示。

圖2 步態(tài)分割算法流程圖Fig.2 Flow chart of gait segmentation algorithm

2.1 步態(tài)序列特征點提取

奶牛一個步態(tài)周期中趾蹄觸地和足底離地時刻前進方向加速度信號均有極值出現，利用小波變換過零點方法提取信號極值點[16-17]，并將該點作為樣本特征點構成樣本序列。信號突變點檢測是對原始信號進行“磨光”的過程，磨光可以理解為對信號進行平滑處理[18-19]。磨光函數計算公式為

(1)

式中θ(t)——低通濾波器

假設θ(t)二次可導，對其求導表示為

(2)

函數ψ(1)(t)和ψ(2)(t)滿足小波容許條件，可作為小波基函數，以此定義的卷積型小波變換為

(3)

式中 *——卷積運算符號s——尺度因子

θs(t)——θ(t)在s下的伸縮

f(t)——原始信號

2.2 模板樣本建立

將連續(xù)步態(tài)序列定義為測試樣本序列，向量表示為Q=(q(1),q(2),…,q(m))，其中m為序列長度。模板樣本向量表示為P=(p(1),p(2),…,p(n))，其中n為序列長度。由于不同跛行程度的奶牛行走模式不同，將步態(tài)樣本分為健康、輕度跛行、中度跛行和嚴重跛行4個組。根據高速攝像機拍攝的視頻手動分割步態(tài)，將每個步幅起點和終點手動標記為模板和黃金標準。為建立能代表大范圍步幅的模板序列，每組手動標記的步幅通過Matlab的interp1函數線性內插至各組樣本中，并取平均值生成模板序列，待分割的測試樣本及生成的模板樣本如圖3所示。

圖3 模板樣本與測試樣本Fig.3 Template sample and test sample

2.3 累積成本矩陣計算

為計算2個不同樣本序列的最小規(guī)整路徑，構造了一個N×M的最小距離矩陣D。該矩陣中的每行表示模板樣本的一個序列到測試樣本的一個序列之間的距離。如果2個樣本序列相似，則局部距離d較小，規(guī)整成本較低。相反，序列不相似，則d較大，規(guī)整成本較大。本文使用歐氏距離計算局部距離，計算公式為

(4)

式中p(n)——模板樣本序列元素

q(m)——測試樣本序列元素

為了尋找模板樣本與測試樣本間最小規(guī)整路徑，定義了一個N×M的累積成本矩陣C。該矩陣表示模板樣本映射到測試樣本各個部分的累積成本，是通過距離矩陣建立累積成本矩陣來搜索最佳規(guī)整路徑的。將距離矩陣D的第1列元素值依次求和得到累積成本矩陣C的第1列元素。而矩陣C的第一行元素由矩陣D的第1行構成，定義該行為起始行，用于計算模板樣本映射到測試樣本的累積成本。累積成本矩陣其余元素計算公式為

C(n,m)=d(n,m)+
min(C(n-1,m-1)，C(n-1,m)，C(n,m-1))

(5)

由式(5)可知，其可以通過將累積成本矩陣C3個方向的最小成本元素與距離矩陣D的各個元素相加來計算。因此，累積成本矩陣中元素表示模板樣本映射到測試樣本的最小累積成本。

2.4 算法改進

由于奶牛前進方向加速度信號極值是隨時間周期性變化的，即每個步態(tài)周期特征點是由時間序列和幅度序列構成的二維序列。本文將該二維序列距離作為改進動態(tài)時間規(guī)整算法的局部距離。設模板樣本的時間為pt，幅值為pa，測試樣本的時間為qt，幅值為qa，因此，樣本序列分別為Ps={(pt(1),pa(1)),(pt(2),pa(2)),…,(pt(n),pa(n))},Qs={(qt(1),qa(1)),(qt(2),qa(2)),…,(qt(m),qa(m))}，將式(4)中的p(n)和q(m)用該二維序列來代替，二維歐氏距離計算公式為

(6)

在將式(6)代入式(5)中，計算模板樣本映射到測試樣本的累積成本矩陣。

2.5 最小規(guī)整路徑確定

最小規(guī)整路徑反映2個樣本相似性的同時可以確定步幅的分割起點及終點。將規(guī)整路徑定義為W，該路徑由成本距離矩陣的每一行元素局部最小值構成，計算公式為

W=(w1,w2,…,wk，…，wK)

(7)

式中K——路徑長度

wk——路徑中的點元素

k——測試樣本序列長度

假設w0為規(guī)整路徑第1行元素，計算式為

w0=min {C(1,1),C(1，2)，…,C(1,m)}

(8)

動態(tài)時間規(guī)整算法對規(guī)整路徑W做了一定限制，只有滿足下面3個條件才可成立。

(1)邊界約束

(9)

任何2個序列都有可能變化，但其各部分的先后次序不會改變，因此所選路徑必定是從網格左下角出發(fā)，在右上角結束。

(2)連續(xù)性約束

對路徑中wk-1=(n′,m′)和wk=(n,m) 2點，需滿足連續(xù)性約束條件

(10)

式(10)表明在計算路徑時不可能跨過某個點去匹配，只能和相鄰點對齊。這可以保證測試樣本和模板樣本中每個坐標都出現在W中。

(3)單調性約束

對路徑中wk-1=(n′,m′)和wk=(n,m) 2點，要求

(11)

式(11)限制了W上面的點必須要隨時間單調進行。由式(10)、(11)可知，對于任意一個網格點(n,m)，下一個通過的網格點只可能是(n,m-1)、(n-1,m-1)和(n-1,m)。最小規(guī)整成本路徑求解過程如圖4所示。

圖4 最小規(guī)整路徑Fig.4 Minimum warping path

通過計算最小規(guī)整路徑從連續(xù)的步態(tài)序列中分割單個步幅。路徑W的起點和終點分別代表步幅的起點和終點，每個分割點的累積成本必須小于分割閾值，因此，選取合適的分割閾值對步態(tài)進行分割。

3 跛行識別模型建立

3.1 步態(tài)特征提取

對于每個分割的步幅分別提取統計學和運動學步態(tài)特征參數，包括均值A[20]、方差[20]、合加速度Ml[20]、運動變化量Mv[20]、信號幅度面積Sma[21]、平均強度Al[21]、偏度Sk[22]、峰度Fk[22]、步態(tài)周期Ts、步幅數Fs、支持時間St和步幅長度Ls。步幅數為連續(xù)序列中分割的步幅數量，支撐時間通過前進加速度數據變化率計算，步態(tài)周期為步幅分割點之間的時間差，而步幅長度由步行總距離除去步幅數得出。

由于健康奶牛和跛行奶牛樣本數量不一致，且可能存在非正態(tài)分布，因此采用Wilcoxon秩和檢驗對提取的特征進行分析，確定跛行奶牛與健康奶牛特征參數間關系。

3.2 跛行識別模型

目前，國內跛行識別研究主要集中在基于計算機視覺和壓力板測量方法上[23-25]，國外主要通過基于各種傳感器及壓敏墊方法采集數據進行跛行識別的研究[26-30]。本文利用邏輯回歸法建立奶牛跛行識別模型[31]，根據統計檢驗結果選取與跛行相關的6個步態(tài)特征參數，建立單因素邏輯回歸模型。將輕度、中度和嚴重度跛行合并為一個類別，稱之為“跛行”。以二分類變量為模型輸出，“1”表示跛行，“0”表示健康。

對于二分類問題可設置一個閾值，如果G(x)≥0.5，則預測y=1，即該奶?；加絮诵?。如果G(x)≤0.5，則預測y=0，即該奶牛為健康。為預測輸出類似于概率一樣介于0～1之間，將sigmoid函數作為映射函數得到邏輯回歸模型,計算公式為

(12)

式中xi——步態(tài)特征值

α——回歸截距β——回歸系數

G——跛行發(fā)生概率

yi——觀測值標簽i——樣本序號

邏輯回歸是非線性模型，利用最大似然估計法對模型系數進行估計。設Gi=G(yi=1|xi)為給定xi條件下得到的跛行yi=1概率，而yi=0的條件概率為G(yi=0|xi)=1-Gi。因此，邏輯回歸模型的似然函數計算公式為

(13)

對式(13)兩邊取對數，得到邏輯回歸對數似然函數，并利用極大似然估計法求最佳回歸系數。

4 結果與分析

4.1 步態(tài)分割結果對比

本文共采集30頭奶牛976個樣本數據，其中健康奶牛(得分1)21頭共613個樣本、輕度跛行(得分2)奶牛2頭共102個樣本、中度跛行(得分3)奶牛4頭共137個樣本、嚴重跛行(得分大于等于4)奶牛3頭共124個樣本。從每組生成模板樣本，計算每組模板樣本與對應測試樣本間累積成本矩陣，根據動態(tài)規(guī)整原則尋找最小規(guī)整路徑，確定分割起點和終點，最終從連續(xù)的步態(tài)序列中分割單個步幅。

圖5為利用前進方向加速度信號分割的單個步幅，步態(tài)周期開始于趾蹄接觸地面終止于趾蹄下一次接觸地面，圖5中紅色虛線之間為一個步態(tài)周期。在趾蹄接觸地面時刻前進方向加速度值反向最大，之后有一個正方向的加速度信號后進入支撐階段，在支撐階段趾蹄相對靜止，因此前進方向加速度值變化不大，圖中黑色虛線之間表示支撐階段。支撐階段后進入擺動階段，擺動階段開始于趾蹄抬起地面，終止于趾蹄接觸地面，趾蹄抬起時刻前進方向加速度值正向最大，圖中黑色虛線和紅色虛線之間為擺動階段。通過分析可知，前進方向加速度信號下峰值均出現在每個步態(tài)周期開始時刻，因此該信號下峰值對應步幅的起點和終點。

圖5 分割的步幅示意圖Fig.5 Schematic of segmented stride

步態(tài)分割主要以減少步幅丟失和最小化錯誤分割為目的。本文利用混淆矩陣來對分割結果進行分析，精確度U、靈敏度V、準確率E計算公式為

(14)

(15)

(16)

式中Tp——正確分割為步幅的樣本數量

Fp——錯誤分割為步幅的樣本數量

FN——無法識別為步幅的樣本數量

TN——無效步幅的樣本數量

對30頭奶牛改進前后的步態(tài)分割結果對比如表1所示。將每個分割的步幅與黃金標準進行比較來驗證分割的準確性。

表1 算法改進前分割結果對比Tab.1 Comparison of segmentation results %

由表1結果可知，本文算法的精確度、靈敏度、準確率平均值可提高到89.53%、95.51%、87.49%，比改進前分別提高了5.31、4.48、8.43個百分點。其中，輕度跛行奶牛步態(tài)分割靈敏度、準確率提高最為明顯，較改進前分別提高7.22、12.75個百分點。而中度跛行奶牛步態(tài)分割精確度提高到86.44%，較改進前提高了10.81個百分點。

為進一步驗證該算法性能，利用峰值檢測法、自相關函數法和本文算法對30頭奶牛進行步態(tài)分割試驗，結果表明，峰值檢測法總準確率87.44%，自相關函數法總準確率88.82%，改進算法在步態(tài)分割準確率方面有所提高，總體準確率提高到90.57%，相較自相關函數法和峰值檢測法分別提高了1.75、3.13個百分點。

4.2 步態(tài)特征提取結果

為研究健康奶牛和跛行奶牛特征參數間關系，從分割的每個步幅提取步態(tài)特征參數，并取均值。利用Wilcoxon秩和檢驗分析各特征值，結果如表2所示，表中特征值用差異性高到低依次進行排列。分析可知，信號幅度面積、平均強度、步態(tài)周期、步幅長度、支撐時間、前進方向加速度均值、運動變化量和方差在健康奶牛和跛行奶牛之間均表現顯著性差異(P<0.05)。而峰度和偏度在2個方向上均不存在顯著性差異。

表2 Wilcoxon秩和檢驗分析結果Tab.2 Analysis results of Wilcoxon rank sum test

4.3 跛行預測模型的測試

根據統計學分析結果，選取支撐時間、步幅長度、平均強度、信號幅度面積、前進方向加速度均值和運動變化量等6個特征，利用邏輯回歸法建立單因素回歸模型。模型輸入為特征值，輸出為二分類結果，其中“1”表示跛行，“0”表示健康，以0.5作為模型分類閾值。將樣本數據的70%作為訓練集來獲得最佳回歸系數，30%作為測試集進行識別模型驗證，通過極大似然估計法分別求出各模型回歸系數，確定邏輯回歸方程。

利用接受者操作特征曲線(Receiver operating characteristic curve，ROC)下的面積[32]評估各預測模型，該面積范圍在0～1之間，值越大表明識別率越高，各預測模型識別率與ROC曲線下的面積見表3。結果可知，單因素邏輯回歸模型中以前進方向加速度均值為自變量的模型獲得較高的識別準確率，為89.45%，ROC曲線下的面積為0.902。其他5個單因素回歸模型識別率均不小于81.72%，ROC曲線下的面積均不小于0.804。

表3 跛行預測模型結果Tab.3 Results of lameness prediction model

5 結論

(1)利用三維加速度傳感器和無線數據采集設備采集奶牛后趾蹄加速度信號，并對其進行預處理，通過分析前進方向加速度信號可以為奶牛步態(tài)的分割提供理論依據。

(2)提出采用改進的動態(tài)時間規(guī)整算法對奶牛步態(tài)進行分割，步態(tài)分割精確度、靈敏度、準確率的均值較改進前分別提高了5.31、4.48、8.43個百分點?？傮w準確率為90.57%，比自相關函數法和峰值檢測法總準確率分別提高了1.75、3.13個百分點。因此，本文算法準確率高，能獲得較好的分割效果。

(3)依據提取的步態(tài)特征，建立跛行識別模型，獲得了較高的識別率，其中以前進方向加速度均值、信號幅度面積、平均強度、運動變化量、支撐時間和步幅長度建立的模型識別率分別為89.45%、86.81%、86.15%、85.71%、83.44%、81.72%。為后期奶牛跛行識別研究奠定了良好的基礎。