基于高斯混合-隱馬爾科夫融合算法識別奶牛步態(tài)時相

摘要：奶牛步態(tài)時相是反映奶牛健康及跛行嚴(yán)重程度的重要指標(biāo)。為準(zhǔn)確自動識別奶牛步態(tài)時相，本研究提出一種融合高斯混合模型 （Gaussian Mixture Model，GMM） 和隱馬爾科夫模型 （Hidden Markov Model， HMM）的無監(jiān)督學(xué)習(xí)奶牛步態(tài)時相識別算法 GMM-HMM。使用慣性測量單元采集奶牛后肢加速度和角速度信號，通過卡爾曼濾波消除噪聲，篩選并提取特征值，構(gòu)建 GMM-HMM模型，實現(xiàn)奶牛靜立相、連續(xù)步態(tài)中的站立相和擺動相等3種步態(tài)時相的自動識別。結(jié)果表明，靜立相識別的準(zhǔn)確率、召回率和 F1分別為89.28%、90.95%和90.91%，連續(xù)步態(tài)中的站立相識別的準(zhǔn)確率、召回率和 F1分別為91.55%、86.71%和89.06%，連續(xù)步態(tài)中的擺動相識別的準(zhǔn)確率、召回率和 F1分別為86.67%、91.51%和89.03%。奶牛步態(tài)分割的準(zhǔn)確率為91.67%，相較于基于事件的峰值檢測法和動態(tài)時間規(guī)整算法準(zhǔn)確率分別提高了4.23%和1.1%。本研究可為下一步基于穿戴式步態(tài)分析的奶牛跛行特征提取提供技術(shù)參考。

關(guān)鍵詞：奶牛跛行;步態(tài)時相;步態(tài)分割;高斯混合模型;隱馬爾科夫模型;卡爾曼濾波

中圖分類號：S24;TP274"""""" 文獻標(biāo)志碼：A"""""""" 文章編號：SA202204003

引用格式：張楷， 韓書慶， 程國棟， 吳賽賽， 劉繼芳.基于高斯混合-隱馬爾科夫融合算法識別奶牛步態(tài)時相[J].智慧農(nóng)業(yè)（中英文）， 2022， 4（2）：53-63.

ZHANG Kai， HAN Shuqing， CHENG Guodong， WU Saisai， LIU Jifang. Gait phase recognition of dairy cows based on Gaussian Mixture model and Hidden Markov model[J]. Smart Agriculture， 2022， 4（2）：53-63.（in Chinese with English abstract）

1引言

奶牛跛行是由于蹄底潰瘍、趾間皮炎和腐蹄病[1] 等各種肢蹄病引起奶牛不適，導(dǎo)致其承重能力不足以及步態(tài)異常。奶牛跛行與飼養(yǎng)方式、牧場環(huán)境衛(wèi)生等有關(guān)[2]，影響奶牛正常生產(chǎn)生活，進而導(dǎo)致產(chǎn)奶量下降、個體產(chǎn)奶年限縮短以及治療成本上升等問題，嚴(yán)重影響牧場經(jīng)濟效益。當(dāng)前中國奶牛肢蹄病發(fā)病率較高，跛行率超過了30%[3]，人工識別存在效率低、主觀判斷存在偏差、勞動力成本較高和容易造成奶牛應(yīng)激反應(yīng)等問題，發(fā)展奶牛跛行自動識別技術(shù)需求強烈[4]。

奶牛跛行步態(tài)分析的關(guān)鍵是在奶牛連續(xù)步態(tài)中識別擺動相和站立相等步態(tài)時相，通過分析奶牛各個步態(tài)時相的運動強度和持續(xù)時間判斷奶牛跛行程度。奶牛的步態(tài)滿足一定的周期性和規(guī)律性，但在現(xiàn)實環(huán)境中，奶牛步態(tài)受環(huán)境和個體差異的影響較大，自適應(yīng)奶牛步態(tài)分析需要針對現(xiàn)實情況采取相應(yīng)的技術(shù)方法。現(xiàn)較多采用計算機視覺技術(shù)和可穿戴式設(shè)備[5]來進行奶牛步態(tài)分析。計算機視覺技術(shù)具有非接觸、免應(yīng)激、低成本、高通量、不干擾奶牛正常活動等優(yōu)點[6]，已被用于奶牛步態(tài)分割[7]、步態(tài)特征提取[8] 以及跛行檢測[9] 中。但計算機視覺容易受到環(huán)境光線、復(fù)雜背景以及奶牛相互遮擋的影響[10]。可穿戴式傳感器[11]具有不限于特定測量區(qū)域以及不受光照和遮擋影響的優(yōu)點，目前，基于可穿戴式傳感器進行奶牛跛行的研究熱點集中在奶牛步態(tài)分割和特征提取方面。常見的奶牛步態(tài)分割方式主要有基于事件的峰值檢測法和基于時域的動態(tài)時間規(guī)整方法兩種。峰值檢測法最先應(yīng)用于奶牛步態(tài)分割和識別[12]，但其容易將一些非步態(tài)的峰值識別為步態(tài)事件，導(dǎo)致低準(zhǔn)確率和高召回率。動態(tài)時間規(guī)整法是指通過計算兩個序列的相似性進行步態(tài)分割，設(shè)定閾值以及建立最優(yōu)匹配模板，自動化程度不高[13]。

隱馬爾科夫模型 （ Hidden Markov Model， HMM ）是馬爾科夫模型發(fā)展而來的機器學(xué)習(xí)算法，核心思想是通過研究觀測序列來辨識不可見狀態(tài)序列。近年來，許多研究在人的步態(tài)識別上使用 HMM 。在人體步態(tài)的研究中分析發(fā)現(xiàn)，步態(tài)信號符合馬爾科夫鏈的時序變化，針對時序數(shù)據(jù)的動態(tài)過程具有魯棒性以及能夠處理牧場條件下奶牛步態(tài)數(shù)據(jù)中的異質(zhì)性，可以進行步態(tài)階段識別[14]。不同步態(tài)分割方法對比顯示，基于 HMM 的分割方法優(yōu)于動態(tài)時間規(guī)整模板匹配法以及峰值檢測法[15]。由于本研究解決的是長時間連續(xù)序列的奶牛步態(tài)時相識別問題，觀測數(shù)據(jù)量較大，所以引入高斯混合模型（ Gaussian Mix ‐ ture Model ，GMM ）。 GMM可以無監(jiān)督地對觀測數(shù)據(jù)進行聚類， GMM 和 HMM 混合算法曾用于人類步態(tài)識別，采用無監(jiān)督的學(xué)習(xí)方式，得到帕金森患者步態(tài)分割，準(zhǔn)確率為96.2%[16]。目前，基于 GMM 和 HMM 融合的奶牛步態(tài)時相識別研究尚未展開。

綜上所述，現(xiàn)有技術(shù)方法初步實現(xiàn)了奶牛步態(tài)分割，但還未深入到奶牛步態(tài)時相識別，并且依舊存在數(shù)據(jù)標(biāo)注困難、準(zhǔn)確率均未達到實際應(yīng)用要求的問題。為進一步改進奶牛步態(tài)時相自動識別方法以及提高其準(zhǔn)確率，本研究提出一種融合 GMM 和 HMM 的無監(jiān)督學(xué)習(xí)奶牛步態(tài)時相識別算法 GMM-HMM ，利用慣性測量單元采集奶牛步態(tài)數(shù)據(jù)，結(jié)合視頻標(biāo)注構(gòu)建奶牛步態(tài)時相識別數(shù)據(jù)集;篩選出冠狀軸角速度的一階差分作為奶牛步態(tài)時相識別的特征值，減少奶牛個體差異對識別效果的影響，并分析奶牛步態(tài)分割效果，為奶牛步態(tài)分析提供一種新的技術(shù)手段。

2材料與方法

2.1試驗設(shè)備與場地

試驗在中國農(nóng)業(yè)機械化科學(xué)研究院北京農(nóng)機試驗站奶牛養(yǎng)殖場進行。試驗對象為5頭隨機選擇的荷斯坦奶牛，利用 MEMS 微型慣性測量單元傳感器（ LPMS-B2型號）采集奶牛后肢加速度和角速度數(shù)據(jù)，采樣頻率為50 Hz ，陀螺儀測量范圍為±2000 dps 。在奶牛活動場內(nèi)，選定一條長20 m 、寬1.5 m走廊作為測量場地，測量場地示意圖如圖1所示。使用佳能 EOS 5D Mark III相機同步拍攝奶牛通過走廊的視頻。

2.2數(shù)據(jù)采集與預(yù)處理

2.2.1數(shù)據(jù)采集

由于奶牛蹄病發(fā)生于后肢的概率較高[17]，以及為方便收集數(shù)據(jù)，統(tǒng)一將傳感器綁在奶牛右后肢上。5頭試驗?zāi)膛Ｅ宕鞣较蛲耆嗤趥鞲衅髦校x x 軸為奶牛前進方向，y 軸為豎直方向;z 軸為冠狀軸方向，如圖2所示。奶牛在測量場地自由活動，相機以50 f/s幀率錄制奶牛通過走廊的視頻。根據(jù)傳感器的時間戳與相機記錄的視頻時間實現(xiàn)圖像數(shù)據(jù)與傳感器數(shù)據(jù)時間同步。試驗中，采用人工引導(dǎo)的方式確保所有試驗?zāi)膛６啻瓮ㄟ^走廊，記錄充足的圖像數(shù)據(jù)。將人工標(biāo)注的視頻數(shù)據(jù)作為奶牛步態(tài)時相識別的黃金標(biāo)準(zhǔn)。試驗采集了2022年2月23日14：00—16：30的奶牛步態(tài)數(shù)據(jù)，每個傳感器有240，000～260，000條數(shù)據(jù)記錄，拍攝的奶牛通過走廊的視頻共計18段。

2.2.2數(shù)據(jù)預(yù)處理

LPMS-B2傳感器數(shù)據(jù)中與奶牛步態(tài)的運動學(xué)數(shù)據(jù)直接相關(guān)的是加速度和角速度，使用3軸加速度和3軸角速度數(shù)據(jù)為基礎(chǔ)進行奶牛步態(tài)時相識別。由于傳感器同時測量9軸數(shù)據(jù)，高頻噪聲較少，噪聲多受系統(tǒng)干擾的影響，為去除系統(tǒng)的噪聲和干擾影響，還原真實的動態(tài)數(shù)據(jù)，采用卡爾曼濾波（Kalman filtering）進行數(shù)據(jù)預(yù)處理。本試驗使用的傳感器集成加速度計、陀螺儀和磁力計三個部分，同時測量 x 、y 、z軸三個方向的線加速度、角速度以及磁場強度。在測量過程中合計9軸的數(shù)據(jù)互相影響并且單位不同，使用卡爾曼濾波可以減輕傳感器噪聲，降低系統(tǒng)誤差和測量誤差。該傳感器內(nèi)置濾波器，在采集數(shù)據(jù)之前已將濾波器設(shè)定完成，得到的數(shù)據(jù)噪聲因素已經(jīng)消除，可以直接利用數(shù)據(jù)進行模型分析。將圖像數(shù)據(jù)采用 FFmpeg軟件進行視頻分幀處理，通過逐幀比對的方式，確定蹄離地和蹄落地兩個步態(tài)事件，再與傳感器數(shù)據(jù)進行匹配。由圖2可知，x軸的加速度可以表示奶牛基本的運動學(xué)特征，z 軸方向為水平且垂直于前進方向向外，z軸角速度可以捕捉奶牛細微的牛蹄運動，可以進一步分析奶牛的步態(tài)情況，所以分析奶牛后肢的x軸加速度和 z軸角速度可以反應(yīng)奶牛步態(tài)事件，更容易實現(xiàn)奶牛步態(tài)事件與傳感器數(shù)據(jù)的同步。x 軸加速度和 z 軸角速度折線圖對應(yīng)的奶牛步態(tài)事件如圖3所示。

2.3訓(xùn)練集和測試集數(shù)據(jù)構(gòu)建

數(shù)據(jù)集構(gòu)建以測試集為基準(zhǔn)，測試集構(gòu)建時使用人工標(biāo)注的方式，將圖像數(shù)據(jù)與傳感器數(shù)據(jù)進行匹配，在傳感器數(shù)據(jù)中標(biāo)出步態(tài)事件（蹄落地，蹄離地），圖像數(shù)據(jù)是隨機采集的，符合測試集構(gòu)建的規(guī)則。觀察6軸傳感器數(shù)據(jù)并做了大量嘗試后發(fā)現(xiàn)，需要找出周期性最明顯的數(shù)據(jù)，方便數(shù)據(jù)標(biāo)注也更容易找到特征值，最終確定將z軸的角速度作為本次試驗數(shù)據(jù)，如圖4所示。

為消除佩戴傳感器時奶牛行為以及剛佩戴好傳感器后奶牛不適應(yīng)帶來的試驗誤差，將原始數(shù)據(jù)的前80，000個記錄點（約26.7 min）的數(shù)據(jù)作為誤差項去除，同理，將后80，000個記錄點同時去除。由于傳感器存在時間戳，利用視頻時間可以實現(xiàn)奶牛步態(tài)時間與傳感器數(shù)據(jù)同步，每頭牛有80～90 min的傳感器數(shù)據(jù)，提取有視頻標(biāo)注的傳感器數(shù)據(jù)作為測試集，剩余的傳感器數(shù)據(jù)去除誤差項作為無監(jiān)督的訓(xùn)練集，其中包含每頭牛至少30 min以上的高質(zhì)量傳感器數(shù)據(jù)，測試集中連續(xù)步態(tài)有162步，訓(xùn)練集中連續(xù)步態(tài)有729步。

2.4 GMM-HMM 融合算法描述

本研究采用 GMM-HMM 融合算法進行奶牛步態(tài)時相識別研究。

奶牛步態(tài)時相為奶牛行走的單個步態(tài)周期中不同階段，本研究將連續(xù)步態(tài)分為擺動相和站立相。擺動相為奶牛蹄離地到蹄著地的過程，站立相為蹄著地到蹄離地的過程。在連續(xù)的步態(tài)序列中，擺動相和站立相是交替出現(xiàn)的，但是本試驗過程中，奶牛除了行走以外，還存在靜止站立的狀態(tài)。因此，本研究將奶牛四肢同時著地靜止不動的狀態(tài)定義為靜立相，區(qū)別于連續(xù)步態(tài)序列中，奶牛四肢交換站立的站立相。靜立相、擺動相和站立相都可以反映奶牛的運動狀態(tài)，同時奶牛步態(tài)時相的識別可以實現(xiàn)對奶牛的步態(tài)分割，識別的連續(xù)站立相和擺動相就是一個完整的步態(tài)。

GMM是多個高斯概率密度函數(shù)的線性組合。為滿足奶牛步態(tài)時相識別的魯棒性，需要奶牛長時間佩戴傳感器，采集的測量數(shù)據(jù)包含靜立相 S1 和連續(xù)步態(tài)中的擺動相 S2和站立相 S3，總共3個步態(tài)時相，每個時相對應(yīng)一個高斯分布，測量數(shù)據(jù)符合這3個時相的高斯混合分布。使用該模型的目的是通過訓(xùn)練找到相應(yīng)參數(shù)θ，進而將高斯混合分布分簇，得到原始的每一個高斯分布，實現(xiàn)類似分類的效果，進而將不同的狀態(tài)分類。高斯混合模型的實現(xiàn)過程即為參數(shù)θ的估計，采用期望最大化（Expectation Maximization ，EM ）算法進行參數(shù)估計[18]，該方法基于極大似然估計算法可表示為：

其中，n是 GMM的階數(shù)，本研究中 n=3;ω i是權(quán)值，且ω i gt;0，ω i =1;Φ i （X ）是第 i 個時相的高斯分布概率密度函數(shù);D 、μ、Σ分別是數(shù)據(jù)的維數(shù)、期望和協(xié)方差。

HMM [19] 由馬爾科夫模型發(fā)展而來，是一種動態(tài)貝葉斯的生成模型以及有向圖模型，由隱狀態(tài)序列和觀測序列構(gòu)成，其中每個觀測數(shù)據(jù)代表一個狀態(tài)。核心思想是由于狀態(tài)序列的不可直接觀測的性質(zhì)，通過研究觀測序列來辨識不可見隱狀態(tài)序列，如圖5所示。

HMM 由三個參數(shù)構(gòu)成：初始狀態(tài)概率π、狀態(tài)轉(zhuǎn)移概率A以及發(fā)射概率（隱狀態(tài)發(fā)生觀測事件的概率） B 。π和A構(gòu)成了隱狀態(tài)之間的狀態(tài)轉(zhuǎn)移，如圖6所示。其中，S1為靜止相，S2、S3為連續(xù)步態(tài)中的擺動相和站立相，aij 為第 i個狀態(tài)轉(zhuǎn)移到第j個狀態(tài)的概率 i，j=1，2，3。B將隱狀態(tài)和觀測序列聯(lián)系，觀測序列由狀態(tài)序列生成。該參數(shù)集可以表示為：

λ={π，A，B }"""""""""""""""" （4）

HMM可以解決三個問題：

（1） 評估問題。給定觀測序列和參數(shù)集λ，計算觀測序列出現(xiàn)的概率，采用前向-后向（ For‐ ward-Backword ）算法實現(xiàn)。

（2） 學(xué)習(xí)問題。給定觀測序列，估計模型參數(shù)集λ，方法的核心是使得觀測序列出現(xiàn)的概率最大的參數(shù)即為模型參數(shù)，采用極大似然估計的方法（ EM算法或者 Baum-Welth 算法）實現(xiàn)。

（3） 解碼問題。已知參數(shù)集λ和觀測序列，計算概率最大的狀態(tài)序列，采用 Viterbi 算法實現(xiàn)。Viterbi算法本質(zhì)上是一個動態(tài)規(guī)劃的最短路問題，找到使得狀態(tài)序列發(fā)射概率最大即由觀測狀態(tài)序列最可能對應(yīng)的狀態(tài)序列。

在本研究中，先進行學(xué)習(xí)問題的訓(xùn)練，利用訓(xùn)練集得出參數(shù)集，再基于得到的模型參數(shù)，利用測試集對觀測序列進行解碼。奶牛的觀測數(shù)據(jù)為觀測序列，奶牛的步態(tài)時相序列為狀態(tài)序列，步態(tài)時相為隱含狀態(tài)。為解決奶牛步態(tài)時相分割的魯棒性，本研究針對的是長時間的連續(xù)觀察序列，采用 GMM 與 HMM 融合模型 （ GMM- HMM ），此時參數(shù)集中的發(fā)射概率 B由一組觀察值的高斯分布密度函數(shù)替代，即發(fā)射概率 B用高斯混合模型的參數(shù)進行替代，新的參數(shù)集為：

λ={π，A，ω，μ，Σ}""""""""""" （5）

其中，π、A分別是 HMM的初始概率和狀態(tài)轉(zhuǎn)移概率;ω、μ和Σ分別為觀測序列高斯混合分布的權(quán)重，均值和協(xié)方差。

將 GMM-HMM 中高斯分布數(shù)設(shè)定為3，對應(yīng)靜立相、擺動相和站立相3種步態(tài)時相。將構(gòu)建好的數(shù)據(jù)集提取特征值，將 z軸角速度的1階差分作為測試集帶入模型，進行學(xué)習(xí)問題。首先對參數(shù)初始概率π和轉(zhuǎn)換概率 A進行初始化，然后根據(jù) EM 算法反復(fù)迭代訓(xùn)練得出參數(shù)ω、μ和Σ，此時已經(jīng)得出3種高斯分布，即 3種狀態(tài)，π直接按照輸入的訓(xùn)練集初始狀態(tài)定義，再根據(jù)應(yīng)用統(tǒng)計的方法計算得出A ，公式如下：

其中，aij 是狀態(tài) i轉(zhuǎn)移到狀態(tài)j的概率; Nij、T 分別是狀態(tài) i轉(zhuǎn)移到狀態(tài)j的個數(shù)、各個狀態(tài)之間轉(zhuǎn)換的總個數(shù)。

GMM-HMM 使用 Python3.8中 hmmlearn.hmm模塊實現(xiàn)。

2.5奶牛步態(tài)時相識別

經(jīng)過對奶牛行走的大量觀察以及查閱資料可知，奶牛的單個步態(tài)周期為1.2 s 左右[20]。為滿足步態(tài)時相識別的自適應(yīng)性和魯棒性，需要對連續(xù)長時間序列的步態(tài)數(shù)據(jù)進行研究。構(gòu)建訓(xùn)練集時，考慮到有監(jiān)督學(xué)習(xí)的訓(xùn)練方式需要人工標(biāo)注的數(shù)據(jù)量過于龐大且效率低，因此本研究采用無監(jiān)督學(xué)習(xí)的訓(xùn)練方式。通過 GMM-HMM 融合算法，從傳感器的運動學(xué)數(shù)據(jù)中找到滿足條件的特征值，對連續(xù)長時間序列數(shù)據(jù)進行無監(jiān)督的聚類以及參數(shù)集訓(xùn)練，得到模型參數(shù)，用測試集驗證結(jié)果，與標(biāo)注數(shù)據(jù)對比，得出識別結(jié)果，確定奶牛步態(tài)時相，并實現(xiàn)奶牛步態(tài)分割。具體流程如圖7所示。

3結(jié)果與分析

3.1特征值選擇

對奶牛步態(tài)的運動學(xué)數(shù)據(jù)進行統(tǒng)計化處理。將2.3中選取的 z 軸角速度采用構(gòu)建滑動時間窗和不使用構(gòu)建滑動時間窗的方式，提取數(shù)據(jù)的最值、均值、方差、標(biāo)準(zhǔn)差、1階差分、線性回歸的斜率和截距，以及平均絕對誤差。最終發(fā)現(xiàn)不加窗的情況下，1 階差分可以很好地區(qū)別奶牛步態(tài)不同時相的數(shù)學(xué)特征。圖8為不同特征值的分布比較情況，截取的是連續(xù)步態(tài)周期的一段，兩側(cè)離散程度大的是擺動相，兩條黑線之間離散程度小的是站立相。由圖8可知，偏度在站立相和擺動相沒有區(qū)分度，最值和均值的重合度過高且沒有區(qū)分兩個相位的直觀特征，標(biāo)準(zhǔn)差區(qū)分度和精度不足，只有1階差分可以直觀且高精度的區(qū)分站立相和擺動相，其中站立相數(shù)據(jù)較為平緩，擺動相數(shù)據(jù)波動較大。

3.2模型參數(shù)訓(xùn)練結(jié)果

經(jīng)過訓(xùn)練得到的模型參數(shù)包括不同時相的初始概率、高斯混合模型的權(quán)重、均值和方差，如表1所示，公式（7） 表示步態(tài)時相的轉(zhuǎn)移概率矩陣A 。轉(zhuǎn)移概率矩陣按照靜立相、站立相、擺動相的順序排列，表示時相之間相互轉(zhuǎn)換的概率。

從所得參數(shù)的結(jié)果可知，初始概率即為訓(xùn)練集的初始時相的概率，模型構(gòu)建時以擺動相為初始輸入時相，滿足隨機性。三個時相的高斯混合模型的權(quán)重都為1，說明每一階段單獨表示相應(yīng)時相，沒有時相之間的交叉影響，提高了模型的可解釋性和準(zhǔn)確性。轉(zhuǎn)移概率的數(shù)值也表明了模型構(gòu)建的準(zhǔn)確性。每一組均值和方差代表一個時相。

3.3奶牛步態(tài)時相識別結(jié)果

參數(shù)訓(xùn)練完成后進入解碼階段。將測試集帶入 GMM-HMM 識別奶牛步態(tài)時相，得到的部分測試結(jié)果如圖9所示。圖中紅色、綠色和藍色部分為奶牛的靜立相、站立相和擺動相。從圖9可以直觀發(fā)現(xiàn)，較為平緩的是奶牛沒有行走的靜立相;站立相是奶牛蹄著地到蹄離地的這一過程，其角速度與圖中體現(xiàn)的一致，從蹄著地開始減少再增加直到蹄離地;連續(xù)步態(tài)的站立相和擺動相交替出現(xiàn)，擺動相的角速度增加到一個步態(tài)周期的最大值，然后減少，直到蹄落地開始受力為止。

采用精確率（ Precision ，P ），召回率（ Re‐ call ， R ），特異率 （Specificity ， Sp），敏感性（Sensitivity ，Se ）和 F1值為評價指標(biāo)。本方法的步態(tài)時相識別結(jié)果中每一條測試集數(shù)據(jù)都會得到一個狀態(tài)，根據(jù)測試視頻標(biāo)注結(jié)果來驗證測試集，結(jié)果如表2所示。

由表2結(jié)果可知，3 種時相在本研究方法的 F1值相差不大，靜立相的 F1值較高為90.91%，識別效果最好，說明靜立相與連續(xù)步態(tài)的兩個時相區(qū)分效果較好;擺動相的召回率和敏感性最高，都為91.51%，說明擺動相的特征較為明顯，識別成功率最高;站立相的準(zhǔn)確率和特異性最高，分別為91.55%和95.51%。3個相位識別結(jié)果的誤差來源主要在于數(shù)據(jù)標(biāo)注時數(shù)據(jù)質(zhì)量以及試驗環(huán)境因素，例如，視頻清晰度不足以準(zhǔn)確判斷奶牛蹄離地、蹄落地這些事件引起的誤差以及搜集數(shù)據(jù)時牛場濕軟地面引起奶牛二次發(fā)力以及站立不穩(wěn)等因素。

通過檢索文獻發(fā)現(xiàn)，大部分針對奶牛步態(tài)的研究關(guān)注奶牛步態(tài)分割方法為提取跛行奶牛步態(tài)時相特征，本研究更進一步研究奶牛步態(tài)時相識別算法。試驗通過人工引導(dǎo)奶牛的方式，采集了奶牛在不同行走速度的步態(tài)慣性測量數(shù)據(jù)，結(jié)果表明本研究提出的方法能夠克服行走速度對步態(tài)時相識別結(jié)果的影響，準(zhǔn)確識別步態(tài)時相，有望實現(xiàn)非約束環(huán)境下奶牛步態(tài)時相的自適應(yīng)識別。

3.4奶牛步態(tài)分割結(jié)果對比分析

針對奶牛連續(xù)步態(tài)階段中的站立相和擺動相進行奶牛步態(tài)分割。在奶牛連續(xù)步態(tài)中，站立相和擺動相構(gòu)成一個完整的步態(tài)周期，根據(jù)連續(xù)步態(tài)中的站立相和擺動相交替出現(xiàn)，奶牛步態(tài)時相識別連續(xù)步態(tài)階段中的站立相和擺動相的結(jié)果即為奶牛步態(tài)分割的結(jié)果。奶牛步態(tài)分割時不再以每條數(shù)據(jù)為標(biāo)準(zhǔn)，使用成功分割的步態(tài)周期為單位進行結(jié)果評價，此時的標(biāo)準(zhǔn)為人工標(biāo)注的完整步態(tài)周期。將成功識別的步態(tài)周期設(shè)定為正樣本，連續(xù)步態(tài)中驗證162步，設(shè)定誤差為5幀（±0.1 s）。為區(qū)別3.3中奶牛步態(tài)時相識別的結(jié)果，重新定義 TP*為準(zhǔn)確識別單個步態(tài)周期的個數(shù)，F(xiàn)P*為非單個步態(tài)周期，但被識別為步態(tài)周期的個數(shù)，F(xiàn)N*為未被成功識別的步態(tài)周期，其中 TP*為121個，F(xiàn)N*為35個，F(xiàn)P*為11個，此時 TN*無意義，得到準(zhǔn)確率為91.67%。

同時，將本研究利用傳感器的奶牛步態(tài)分割方式準(zhǔn)確率結(jié)果與基于圖像處理的奶牛步態(tài)分割算法，包括支持向量機、K最鄰近法以及隨機森林等進行對比。峰值檢測法奶牛步態(tài)分割準(zhǔn)確率為87.44%，動態(tài)時間規(guī)整奶牛步態(tài)分割準(zhǔn)確率為90.57%[13]，本研究奶牛步態(tài)分割準(zhǔn)確率為91.67%，不僅識別了奶牛步態(tài)時相，而且相較于峰值檢測法和動態(tài)時間規(guī)整分別將奶牛步態(tài)分割準(zhǔn)確提高4.23%和1.1%。與基于圖像處理的奶牛步態(tài)分割算法支持向量機、K最近鄰和隨機森林的準(zhǔn)確率分別為87.60%、83.90%和82.56%，本研究方法分別提高了4.07%、7.7%和9.11%。

本試驗步態(tài)分割成功識別160步，蹄落地的總誤差為-14.9 s ，平均每步誤差0.093 s ，蹄離地的總誤差為-7.2 s ，平均每步誤差-0.045 s ，滿足0.1 s的誤差。其中“-“表示滯后（視頻和傳感器都為50 Hz）。設(shè)立的正樣本是站立相，站立相是奶牛蹄著地到奶牛蹄離地這一過程，蹄著地時奶牛會出現(xiàn)站立不穩(wěn)、滑行，以及二次用力重新站穩(wěn)的現(xiàn)象，導(dǎo)致蹄著地這一過程誤差較大（即站立相開始階段），模型對蹄著地這一事件判別不準(zhǔn)確，蹄落地后，傳感器會再次收到抖動的信號，直到穩(wěn)定狀態(tài)，才會識別站立相，進而導(dǎo)致蹄著地這一事件識別會滯后。雖然蹄離地這一事件，不會有二次受力的情況，誤差會大幅減小，但是松軟、粘連的地面情況，也會出現(xiàn)較小的滯后情況。

4結(jié)論

本研究利用可穿戴式慣性測量單元采集奶牛后肢的3軸加速度和3軸角速度數(shù)據(jù)，經(jīng)過數(shù)據(jù)預(yù)處理和特征值甄別，選擇 z軸角速度的1階差分作為特征值，通過 GMM-HMM融合算法進行奶牛步態(tài)時相的識別。結(jié)果為奶牛步態(tài)時相和奶牛步態(tài)分割提供了理論依據(jù)，具體結(jié)論如下。

（1） 提出了一種 GMM-HMM 融合的無監(jiān)督學(xué)習(xí)奶牛步態(tài)時相識別算法，識別出了奶牛的靜立相、連續(xù)步態(tài)階段的擺動相和站立相，其中靜立相識別的準(zhǔn)確率、召回率和 F1分別為89.28%、90.95%和90.91%，連續(xù)步態(tài)中的站立相識別的準(zhǔn)確率、召回率和 F1分別為91.55%、86.71%和89.06%，連續(xù)步態(tài)中的擺動相識別的準(zhǔn)確率、召回率和 F1分別為86.67%、91.51%和89.03%。進一步識別了奶牛步態(tài)時相，且識別效果較好，為后期奶牛跛行特征提取提供了技術(shù)方法。

（2） 根據(jù)連續(xù)步態(tài)的站立相和擺動相實現(xiàn)了奶牛步態(tài)分割，準(zhǔn)確率達到了91.67%，相較于基于事件的峰值檢測法和動態(tài)時間規(guī)整算法準(zhǔn)確率分別提高了4.23%和1.1%，分割效果較好。

為滿足非約束條件下步態(tài)識別的應(yīng)用要求，今后將重點測試奶牛在不同的行走速度、跛行狀態(tài)以及地面類型的條件下的模型性能，以提升模型的魯棒性，并結(jié)合跛行引起的奶牛躺臥、采食等日常行為波動狀況分析，提高跛行判別的準(zhǔn)確性。

參考文獻：

[1] TADICH" N，" FLOR" E，" GREEN" L. Associations" be‐tween hoof lesions and locomotion score in 1098 un‐ sound dairy cows[J]. The Veterinary Journal， 2010， 184（1）：60-65.

[2] OEHM A W，" KNUBBEN-SCHWEIZER" G，" RIEGERA， et al. A systematic review and meta-analyses of risk factors" associated" with" lameness" in" dairy" cows[J]. BMC Veterinary Research， 2019， 15（1）： ID 346.

[3] 邵軍， 林為民， 孫新文， 等.石河子地區(qū)奶牛肢蹄病類型和發(fā)病情況調(diào)查分析[J].中國奶牛， 2017（1）：30-33.

SHAO J， LIN W， SUN X， et al. Investigation and anal‐ ysis of limb hoof disease types and incidence in dairy cows in Shihezi area[J]. China Dairy Cattle， 2017（1）：30-33.

[4] 韓書慶， 張晶， 程國棟， 等.奶牛跛行自動識別技術(shù)研究現(xiàn)狀與挑戰(zhàn)[J].智慧農(nóng)業(yè)（中英文）， 2020， 2（3）：21-36.

HAN S， ZHANG J， CHENG G， et al. Research status and" challenges" of automatic" identification technology for cow lameness[J]. Smart Agriculture， 2020， 2（3）：21-36.

[5] 程國棟， 吳建寨， 邢麗瑋， 等.基于 IMU的細粒度奶牛行為判別[J].中國農(nóng)業(yè)大學(xué)學(xué)報 ， 2022， 27（4）：179-186.

CHENG G， WU J， XING L， etc. The behavior discrimi‐ nation" of fine-grained" dairy" cows" based" on" IMU[J]. Journal of China Agricultural University， 2022， 27（4）：179-186.

[6] LI Z， ZHANG Q， LYU S， et al. Fusion of RGB， opticalflow" and" skeleton" features" for the" detection" of lame‐ ness in dairy cows[J]. Biosystems Engineering， 2022， 218：62-77.

[7] ALSAAOD M， NIEDERHAUSER J J， BEER G， et al.Development and validation of a novel pedometer algo‐ rithm to quantify extended characteristics of the loco‐ motor behavior of dairy cows[J]. Journal of Dairy Sci‐ ence， 2015， 98（9）：6236-6242.

[8] ZHAO K， BEWLEY J M， HE D， et al. Automatic lame‐ness detection in dairy cattle based on leg swing analy‐ sis with an image processing technique[J]. Computers and Electronics in Agriculture， 2018， 148：226-236.

[9] 宋懷波， 姜波， 吳倩， 等.基于頭頸部輪廓擬合直線斜率特征的奶牛跛行檢測方法[J].農(nóng)業(yè)工程學(xué)報 ，2018， 34（15）：190-199.

SONG H， JIANG B， WU Q， et al. Detection method ofcow lameness based on fitting straight line slope char‐acteristics of head and neck contour[J]. Transactions ofthe CASE， 2018， 34（15）：190-199.

[10] JIANG" B，" SONG" H，" HE" D. Lameness" detection" ofdairy cows based on a double normal background sta‐tistical model[J]. Computers" and" Electronics" in Agri‐culture， 2019， 158：140-149.

[11] PENG Y， KONDO N， FUJIURA T， et al. Classificationof multiple" cattle behavior patterns using" a recurrentneural network with long short-term memory and iner‐tial measurement units[J]. Computers" and Electronicsin Agriculture， 2019， 157：247-253.

[12] HALADJIAN J， BRüGGE B， NüSKE S. An approachfor" early" lameness" detection" in" dairy" cattle[C]// The2017 ACM" International" Joint" Conference" on" Perva‐sive" and Ubiquitous" Computing" and" The 2017 ACMInternational" Symposium" on" Wearable" Computers.New York， USA： Association for Computing Machin‐ery， 2017.

[13]蘇力德， 張永， 王健， 等.基于改進動態(tài)時間規(guī)整算法的奶牛步態(tài)分割方法[J].農(nóng)業(yè)機械學(xué)報， 2020， 51（7）：52-59.

SU L， ZHANG Y， WANG J， et al. Gait segmentationmethod" for" dairy" cows" based" on "improved" dynamictime warping algorithm[J]. Transactions of the CSAM，2020， 51（7）：52-59.

[14]張向剛， 唐海， 付常君， 等.一種基于隱馬爾科夫模型的步態(tài)識別算法[J].計算機科學(xué) ， 2016， 43（7）：285-289， 302.

ZHANG X， TANG H， FU C， et al. A gait recognitionalgorithm based on hidden Markov model[J]. Comput‐er Science， 2016， 43（7）：285-289， 302.

[15] HAJI" GHASSEMI N，" HANNINK" J，" MARTINDALEC F， et al. Segmentation of gait sequences in" sensor-based movement analysis： A comparison of methods inParkinson's" disease[J]. Sensors （Basel）， 2018， 18（1）：ID 145

[16] ROTH N， KUDERLE A， ULLRICH M， et al. HiddenMarkov model based stride segmentation on unsuper‐vised" free-living" gait" data" in" Parkinson's" disease" pa‐tients[J]. Journal of Neuro Engineering and Rehabilita‐tion， 2021， 18（1）： ID 93.

[17]王東源， 王玉舜， 谷禹， 等.天津地區(qū)奶牛蹄病分析研究[J].天津農(nóng)學(xué)院學(xué)報， 2015， 22（2）：37-39.

WANG" D，" WANG Y，" GU Y，" et" al. Analysis" of cowhoof disease in Tianjin[J]. Journal of Tianjin Agricul‐tural University， 2015， 22（2）：37-39.

[18] XIA" R，" ZHANG" Q，" DENG" X. Multiscale" Gaussianconvolution" algorithm" for" estimate" of Gaussian" mix‐ture" model[J]. Communications" in" Statistics-Theoryand Methods， 2018， 48（23）：5889-5910.

[19] RABINER L R. A Tutorial on Hidden Markov modelsand" selected" applications" in" speech" recognition[M]// WAIBEL A， LEE K-F. Readings in speech recognition. San Francisco： Morgan Kaufmann， 1990.

[20] ALSAAOD M， LUTERNAUER M， HAUSEGGER T，et al. The cow pedogram—Analysis of gait cycle vari‐ables allows the detection of lameness and foot pathol‐ogies[J]. Journal of Dairy Science， 2017， 100（2）：1417-1426.

Gait Phase Recognition of Dairy Cows based on Gaussian Mixture Model and Hidden Markov Model

ZHANG Kai， HAN Shuqing， CHENG Guodong， WU Saisai， LIU Jifang*

（Agricultural Information Institute of Chinese Academy of Agricultural Sciences/ Key Laboratory of AgriculturalBlockchain Application， Ministry of Agriculture and Rural Affairs， Beijing 100081， China ）

Abstract： The gait phase of dairy cows is an important indicator to reflect the severity of lameness. IThe accuracy of available gait segmentation methods was not enough for lameness detection. In this study， a gait phase recognition method based on Gaussian mixture model （GMM） and hidden Markov model （HMM） was proposed and tested. Firstly， wearable inertial sensors LPMS-B2 were used to collect the acceleration and angular velocity signals of cow hind limbs. In order to remove the noise of the system and restore the real dynamic data， Kalman filter was used for data preprocessing. The first-order difference of the an‐ gular velocity of the coronal axis was selected as the eigenvalue. Secondly， to analyze the long-term continuous recorded gait se‐ quences of dairy cows， the processed data was clustered by GMM in the unsupervised way. The clustering results were taken as the input of the HMM， and the gait phase recognition of dairy cows was realized by decoding the observed data. Finally， the cow gait was segmented into 3 phases， including the stationary phase， standing phase and swing phase. At the same time， gait segmentation was achieved according to the standing phase and swing phase. The accuracy， recall rate and F1 of the stationary phase were 89.28%， 90.95% and 90.91%， respectively. The accuracy， recall rate and F1 of the standing phase recognition in con‐ tinuous gait were 91.55%， 86.71% and 89.06%， respectively. The accuracy， recall rate and F1 of the swing phase recognition in continuous gait were 86.67%， 91.51% and 89.03%， respectively. The accuracy of cow gait segmentation was 91.67%， which was 4.23% and 1.1% higher than that of the event-based peak detection method and dynamic time warping algorithm， respec‐ tively. The experimental results showed that the proposed method could overcome the influence of the cow's walking speed on gait phase recognition results， and recognize the gait phase accurately. This experiment provides a new method for the adaptive recognition of the cow gait phase in unconstrained environments. The degree of lameness of dairy cows can be judged by the gait features.

Key words： dairy cow lameness; gait phase; gait segmentation; gaussian mixture model; hidden Markov model; Kalman filtering