基于聲學沖激響應的奶山羊呼吸及反芻同步監測方法

王天本 劉現濤 李張本 嚴宏昊 宋懷波 胡 瑾

(1.西北農林科技大學機械與電子工程學院, 陜西楊凌 712100;2.農業農村部農業物聯網重點實驗室, 陜西楊凌 712100)

0 引言

近年來羊奶制品的市場規模不斷擴大。據統計,2020年全國羊奶市場規模為160億元,2025年將突破200億元[1]。隨著羊奶市場的蓬勃發展,智能化是奶山羊產業發展的必然趨勢。奶山羊個體生理體征監測是智能化養殖不可或缺的環節[2]。呼吸及反芻是奶山羊基本的生理體征,準確地獲取奶山羊個體的呼吸和反芻信息對于評估奶山羊健康狀況,及時發現奶山羊患病個體具有重要意義[3]。

傳統的呼吸和反芻監測以人工觀測為主。人工觀察并計量胸脯起伏、咀嚼次數[4-5],效率低下且無法長時間持續監測。隨著信息技術和傳感技術的迅速發展,自動化呼吸及反芻監測成為奶山羊智能化養殖領域的研究熱點。現有的自動化奶山羊呼吸及反芻監測技術可分為接觸式和非接觸式兩大類。接觸式呼吸監測方法利用穿戴式設備(如壓敏、氣敏或熱敏傳感器等)測定胸脯部運動、呼吸聲、呼吸氣流等相關參數實現呼吸監測[6-8];接觸式反芻監測方法主要以加速度傳感器、壓力傳感器測量下顎移動狀態實現反芻監測[9-12]。雖然上述方法可以實現呼吸和反芻的可靠監測,但需要在牲畜身體特定部位穿戴或植入傳感器,干擾牲畜自然狀態。此外,穿戴式設備易被刮蹭、啃咬、尿液浸泡而出現損壞,實施和維護成本較高。

現有的非接觸式牲畜呼吸和反芻監測方法主要以視頻圖像[13-18]、激光測距[19]、熱成像[20]、聲音識別[21-26]等技術為基礎構建。基于視頻圖像的方法雖能進行準確監測,但受可視角度限制,對被測牲畜的朝向敏感,當被測牲畜背對設備時系統無法監測呼吸和反芻。此外,基于視頻的方法易受環境光線影響干擾。基于熱成像的監測方法雖不受光線影響,但設備價格昂貴,難以推廣應用。基于激光測距的方法由于受激光信號強指向性的影響對被測牲畜的朝向敏感。現有基于聲音識別技術的方法,通過拾音器近距離采集牲畜呼吸過程中鼻腔發出的聲音和反芻時上下頜撞擊產生的聲音,利用機器學習算法實現呼吸和反芻監測。此類方法依賴大量離線數據,對環境噪聲敏感且距離分辨率不足。此外,上述幾類方法都無法實現牲畜呼吸和反芻的同步監測。近年來研究人員嘗試將超聲用于人的非接觸式呼吸監測[27-28],并獲得了良好的呼吸監測效果。與其它監測技術相比,基于超聲的監測技術范圍廣、不受光線的影響,為奶山羊的呼吸及反芻監測提供了理論依據。

沖激響應[29-30]作為衡量系統響應隨時間變化的工具,可以描述聲學信號到達拾音器的延遲及衰減情況,具備良好的時間分辨特性。鑒于上述分析,本文提出一種基于聲學沖激響應的靜臥奶山羊呼吸及反芻同步監測方法。通過聲學沖激響應捕捉呼吸和反芻引起的多徑信號的變化,實現對單只靜臥奶山羊朝向魯棒的呼吸及反芻同步監測,以期為奶山羊患病個體及其它特殊個體的生理體征監測提供技術支撐。

1 聲學信號捕獲與數據采集

1.1 室內聲學回波模型

聲波信號本質上為機械波,無法穿透墻體和室內的其它物體,因此在室內存在豐富的聲學多徑反射效應。在圈舍養殖場景中布設揚聲器和拾音器,揚聲器與拾音器以類似雷達的方式同步工作。揚聲器發出預設的超聲信號,經過室內的一次或多次反射,最終被拾音器接收。將圈舍環境及被測羊只看作一個系統,利用收發信號估計系統聲學沖激響應,刻畫反射信號幅度及時延隨時間的變化。奶山羊呼吸過程中的胸脯起伏和反芻過程中的嘴部咀嚼動作均為周期性運動,導致經過其反射的聲波信號的幅度和延遲發生周期性變化。然而,拾音器接收到的信號不僅包含奶山羊呼吸及反芻運動的反射信號,還包含室內的其它反射信號,為了從眾多反射信號中提取奶山羊呼吸及反芻運動的反射信號,本文構建了室內回波模型。如圖1所示,拾音器接收的回波主要包括經過奶山羊胸脯部的反射信號(圖中綠色箭頭),經過奶山羊嘴部反射的信號(圖中黃色箭頭),經過奶山羊胸脯部及嘴部的反射信號(圖中紅色箭頭),羊的其它部位反射信號(圖中藍色箭頭),靜態環境反射信號(圖中紫色箭頭),環境噪聲(圖中白色箭頭)。

圖1 室內聲學回波模型

綜上,在室內環境中,拾音器所接收到的回波xr(n)可以表示為

xr(n)=xbre(n)+xchew(n)+xmult(n)+xdyn(n)+

xsta(n)+N(n)

(1)

式中xbre(n)——經過奶山羊胸脯部反射,最終被拾音器接收的信號,包括直接經過奶山羊胸脯部反射信號xbre_dir(n),以及同時經過周圍環境(如墻壁、門、窗等)和奶山羊胸脯部的反射信號xbre_indir(n)

xchew(n)——經過奶山羊嘴部反射,最終被拾音器接收的信號,包括直接經過奶山羊嘴部反射信號xchew_dir(n)以及同時經過周圍環境(如墻壁、門、窗等)和奶山羊嘴部的反射信號xchew_indir(n)

xmult(n)——經過奶山羊胸脯部及嘴部的反射信號,包括直接反射信號xmult_dir(n)和間接反射信號xmult_indir(n)

xdyn(n)——經過奶山羊其它部位(如耳部、腿部、尾部等)一次或多次反射,最終被拾音器接收到的回波信號

xsta(n)——經過靜態環境的一次或多次反射信號

N(n)——圈舍環境中的噪聲

奶山羊在圈舍內處于靜息狀態的位置和朝向不固定,當奶山羊背對設備或被遮擋時,xbre(n)、xchew(n)和xmult(n)中不存在或存在很少的直接反射信號,但由于室內環境中存在豐富的聲學多徑反射,因此xbre_indir(n)、xchew_indir(n)和xmult_dir(n)一直存在。在實際圈舍環境中,由于奶山羊胸脯部的面積遠小于靜態環境的反射面積,因此xbre(n)、xchew(n)和xmult(n)中包含反射路徑的數量和信號能量遠小于xsta(n)中包含反射路徑的數量和信號能量。但由于奶山羊呼吸過程中的胸脯起伏和反芻過程中的嘴部咀嚼動作均為周期性運動,靜態環境反射信號xsta(n)沒有經過奶山羊的胸脯部及嘴部反射,不包含奶山羊的運動信息,因此可以通過周期性分析反射信號來提取xbre(n)、xchew(n)和xmult(n)。提取上述周期性變化的反射信號,利用呼吸頻率和反芻頻率的差異[31-32]實現呼吸信號和反芻信號的分離。

1.2 回波音頻采集

于2022年6—7月在西北農林科技大學畜牧教學試驗基地隔離圈舍(3 m×2.5 m×2.9 m)進行靜臥奶山羊回波音頻采集,選取具備反芻能力的7只健康西農薩能奶山羊,分別進行單獨試驗。為了避免羊只對環境改變產生的生理變化,每只奶山羊都提前3 d在隔離圈舍內作環境適應性飼養。試驗布置如圖2所示,在圈舍內布設聲波多路同步收發平臺,聲波收發平臺由4個寬頻揚聲器和嵌有兩個全向拾音器的聲卡組成。聲波收發平臺通過USB接口與便攜式計算機相連,呼吸及反芻監測算法在Matlab軟件上運行。為了音頻采集的效果及避免被羊只破壞,聲波收發平臺放置在圈舍左側,在外放置保護柵欄。試驗采集數據以音頻幀的形式存儲到硬盤中,數據格式為.mat,幀率為20 f/s。每段音頻數據時長為20 min。

圖2 音頻采集試驗布署

試驗平臺為一臺Windows 10專業版系統的便攜式計算機,處理器配置為Intel(R)Core(TM)i7-4720 HQ,主頻為2.60 GHz,內存為8 GB,算法開發平臺版本為Matlab 2020a。收發平臺采樣頻率fs為96 kHz,以頻率變化連續的正弦調頻連續波作為發射信號,參數設置為:載波頻率fc為38 kHz,掃頻帶寬B為2 kHz,掃頻周期T為0.02 s,一個信號周期T′為0.05 s。4個揚聲器的總功率為50 W,揚聲器靈敏度為94 dB。當fc=38 kHz時,揚聲器的波束寬度為±14°,拾音器的靈敏度、信噪比和總諧波畸變分別為-26 dBFS、64.3 dB和0.2%。

1.3 音頻選擇

對采集到的圈舍內靜臥奶山羊音頻數據進行篩選,去除掉奶山羊進食、未發生反芻行為、奶山羊活動等無效音頻數據后,共得到7段奶山羊試驗音頻。如表1所示,試驗中奶山羊相對聲波收發平臺朝向涉及正朝、側朝和背朝3種方向,且存在雨聲、其它羊叫聲、音樂聲等噪聲。7段測試數據的時長、相對聲波收發平臺朝向均具有較大差異,且存在其它干擾因素,可用于驗證本文方法的適用性。

表1 靜臥奶山羊音頻信息

2 奶山羊呼吸及反芻監測算法

算法主要由3部分組成:沖激響應序列估計,用于定量描述每一路多徑信號的變化;呼吸及反芻信號分離,通過自相關函數度量沖激響應時序序列,根據呼吸頻率和反芻頻率的差異分離呼吸和反芻信號;呼吸及反芻波形可視化,將分離得到的呼吸和反芻信號,經幅度歸一化、相位同步后,分別合成呼吸波形和反芻波形。算法整體方案如圖3所示。

圖3 算法整體方案

2.1 沖激響應序列估計

為了準確地估計圈舍環境中的系統沖激響應,確保沖激響應估計的準確性并避免對羊只造成干擾,需要設計合適的發送信號。本文采用載波頻率大于37 kHz(羊聽覺范圍上限)的連續調頻波作為發送信號。由于常見的鋸齒調頻或者三角調頻均存在頻率突變點,發送信號存在頻譜泄漏現象,導致揚聲器發出的聲波混雜有羊聽覺范圍內的噪聲。本文選用頻率變化連續的正弦調頻構建發送信號。其瞬時頻率可表示為

(2)

發送信號在時刻t的瞬時相位u(t)可表示為

(0≤t≤T)

(3)

發送信號表示為

xs(t)=cos(u(t))=

(0≤t≤T)

(4)

為實現持續的呼吸及反芻監測,需要連續發送xs(t),拾音器同步接收當前時刻反射回波。但直接發送xs(t)會導致回波信號存在混疊現象,即某一時刻拾音器接收到的數據為多次發送信號回波的疊加,不滿足系統沖激響應估計的零狀態要求,導致沖激響應的估計錯誤。為了使每次拾音器接收到的數據為一次發送信號的完整回波,需要為發送信號加入一定的占空比,公式為

(5)

式中w(t)——窗函數,用于緩解信號幅度突變導致的頻譜泄漏

最終發送信號的離散形式為

xs(n)=

(6)

式中Ts——采樣間隔

h(n)——圈舍聲學沖激響應

為了避免接收回波出現混疊,T′-T應大于圈舍最大多徑時延。通常圈舍空間越大,最大的多徑時延也就越大。圖4a、4b分別為發送信號的時域波形和時頻域波形(fs=96 kHz,fc=38 kHz,B=2 kHz,T=0.02 s,T′=0.05 s)。

圖4 發送信號時域與時頻域圖

為了減少環境噪聲干擾,收到的回波濾除發送掃頻頻段周圍500 Hz頻段的信號,以消除發送信號頻帶以外的噪聲干擾。去除噪聲后,拾音器接收到的回波模型為

xr(n)=h(n)?xs(n)

(7)

根據卷積定理可得

(8)

式中H(n)——h(n)的傅里葉變換

Xr(n)——xr(n)的傅里葉變換

Xs(n)——xs(n)的傅里葉變換

(9)

ε(n)——正則項,一般設計為在xs(n)的頻段內取很小的值而在xs(n)頻段外取較大的值以屏蔽環境噪聲

最終沖激響應可估計為

(n)=IFFT(Xs(n)C(n))

(10)

以圖4中所示例子為發送信號,圖5為某圈舍內環境系統沖激響應的估計結果。橫軸表示每個多徑反射信號的時間延遲,縱軸表示多徑信號幅度。

圖5 系統沖激響應序列時域圖

[S1S2…Sk]

(11)

式中Si——時間延遲為i/fs的多徑信號的幅度變化序列

M——分析窗口長度

由于聲波收發設備不是嚴格的線性系統,Si會隨著時間的推移產生非線性趨勢,因而要對Si進行去趨勢處理。

2.2 呼吸及反芻信號分離

呼吸過程中的胸脯起伏和反芻過程中的嘴部咀嚼動作均具有周期性,經過胸脯部及嘴部的直接或間接多徑反射信號,均會出現周期性波動。因此,通過度量Si的周期性即可實現直接和間接呼吸及反芻動作反射信號的提取。

自相關函數是進行周期性度量的良好工具,Si的自相關函數定義為

(12)

(13)

式中cn——Si的自協方差

圖6a給出了4個不同的Si,分別表示:僅經過胸部和靜態環境反射的多徑信號(黑色);僅經過嘴部和靜態環境反射的多徑信號(紅色);同時經過嘴部、胸部及靜態環境反射的多徑信號(藍色);未經過胸部及嘴部反射的多徑信號(黃色)。可以看到僅經過胸部、僅經過嘴部和同時經過嘴部和胸脯反射的多徑信號均具有一定的周期性,而未經過嘴部和胸部反射的信號無周期性。圖6b為上述4個Si的自相關函數值。可以看到,有周期性的Si的自相關函數呈現出多個明顯的峰值且峰值之間間隔均勻。假設某一時刻Si序列的自相關函數前i個峰值的幅度為peaki,對應的峰值索引位置為li,0≤i≤m,令自相關函數峰值間隔linr=li+1-li,因此,基于自相關函數的挑選周期性Si的規則可定義為

圖6 序列與自相關函數值

peakmax>PeakheightThrd

其中,peakmax表示自相關函數峰值的最大值;PeakheightThrd表示自相關函數峰值幅值的閾值;PeriodThrd表示自相關函數周期的閾值;MidD表示奶山羊完成一次反芻的平均時間。

基于上述規則即可提取出所有包含呼吸和反芻信息的Si。

基于自相關函數提取出來的周期性Si可能同時包含呼吸和反芻信息(如圖6a中藍曲線所示)。由于奶山羊呼吸頻率和反芻頻率差異明顯,因此可以通過構建帶通濾波器組分離呼吸和反芻信號。為了適應不同羊只的呼吸和反芻頻率,首先計算所有周期性Si的幅度譜,如圖7所示,根據幅度譜峰值確定帶通濾波器的截止頻率,然后使用構造的帶通濾波器進行濾波。圖8展示了圖6a中同時包含呼吸和反芻信息的Si經過帶通濾波器組分離得到的呼吸波和反芻波。

圖8 濾波前后序列

2.3 呼吸及反芻運動波形可視化

從多個不同的周期性Si可以分離出多路呼吸信號和多路反芻信號。融合多路呼吸信號(反芻信號)即可得到呼吸波(反芻波)。由于不同的多徑信號有不同的幅度和反射路徑,從不同的周期性Si分離得到的多個呼吸信號(反芻信號)幅度不同且可能出現相位差為π的現象。圖9a和圖9e分別為從不同周期性Si分離得到的多個呼吸信號和多個反芻信號。直接將這些呼吸信號(反芻信號)疊加會使幅度較小的信號被淹沒,相位相反的信號相互抵消。為解決上述問題,首先將所有呼吸信號(反芻信號)進行幅度歸一化,然后對相位進行同步,最后進行疊加平均(詳見算法1)。圖9b和圖9f分別為歸一化后的多個呼吸信號和多個反芻信號,圖9c和圖9g分別為相位同步后的多個呼吸信號和多個反芻信號,圖9d和圖9h為最終合成的該時刻的呼吸信號和反芻信號。

圖9 波形合成

算法1:呼吸和反芻波形合成算法

輸出:當前時刻合成的呼吸波Wbcurr及反芻波Wccurr。

Bwcurr←Synthesizer(Sb,Wblast);

Cwcurr←Synthesizer(Sc,Wclast);

Procedure Synthesizer(S,W)

forifrom 2 toNdo

else

end if

end for

Wcurr←Wcurr./N∥對合成的波形的幅度進行求和平均處理

if ‖Wlast+Wcurr‖1<‖Wlast-Wcurr)‖1do∥ 與前一時刻相位校正

Wcurr←-Wcurr

end if

算法首先將當前呼吸(反芻)波初始化為歸一化后的第一個呼吸(反芻)信號。然后逐個檢查其它呼吸(反芻)信號,若與當前呼吸(反芻)波相位相近則將其直接疊加到當前呼吸(反芻)波,否則沿時間軸反轉后疊加到當前呼吸(反芻)波。然后對當前呼吸(反芻)波做平均處理。最后參考上一時刻的呼吸(反芻)波對當前呼吸(反芻)波相位進行校正,確保前后兩次合成的呼吸波相位連續。

3 試驗結果與分析

3.1 不同個體奶山羊呼吸及反芻監測結果分析

表1中7段音頻數據采集時奶山羊的位置、對應的某一幀呼吸及反芻波形如圖10所示。可以看到即使被測羊只背對著聲波收發設備,系統仍能監測到奶山羊的呼吸及反芻運動。

采用本文呼吸及反芻監測方法對7段音頻數據進行驗證,結果如表2所示。表2中呼吸和反芻次數的真實值由人工觀察得到;呼吸和反芻次數的計算值由監測結果的波形數得到。通過定義相對誤差,可以評估系統計算值的可信度。漏幀率為漏幀數與總幀數的比值。由表2可知,由監測算法計算的呼吸及反芻次數與真實值相比,具有較高的一致性,其中呼吸平均相對誤差為2.60%,呼吸最大相對誤差為6.67%,反芻平均相對誤差為3.51%,反芻最大相對誤差為4.76%,平均漏幀率為2.49%。綜上結果表明,本文方法能同步監測奶山羊的呼吸及反芻行為,對靜臥奶山羊具有較高的跟蹤精度。

表2 靜臥奶山羊呼吸次數、反芻次數、相對誤差和漏幀率

從表2中可以看出,不同奶山羊的呼吸頻率和反芻頻率具有差異。以表1中第Ⅵ、Ⅴ段的音頻數據(對應圖10第Ⅵ、Ⅴ號奶山羊)為例,呼吸及反芻追蹤結果繪制的時間-頻率曲線如圖11所示。通過觀察該曲線可以發現,即使在奶山羊呼吸頻率及反芻頻率相差較大時,系統仍能穩定追蹤。

圖11 Ⅴ、Ⅵ號奶山羊時間-頻率曲線

3.2 奶山羊朝向對監測結果的影響

奶山羊不同朝向下的統計結果如表3所示,從表3中可以看出,奶山羊側朝聲波收發設備與正朝設備的監測結果沒有明顯差別。Ⅲ、Ⅵ號羊背對聲波收發設備時漏幀率略有上升,平均漏幀率為5.05%。但相較于其它朝向,Ⅲ、Ⅵ號羊呼吸平均相對誤差為3.34%,反芻平均相對誤差為3.21%,沒有明顯的呼吸和反芻相對誤差偏差,滿足實際圈舍環境下監測的需求。這是因為室內環境中存在豐富的聲學多徑反射效應,如1.1節中所述,即使被測羊只背朝聲波收發設備,拾音器也總能接收到經過被測羊只胸脯部和嘴部間接反射的聲波信號,因此系統對被測羊只的朝向具有較強的魯棒性。

表3 靜臥奶山羊不同朝向下的相對誤差和漏幀率

漏幀率略有上升的推斷原因是,當奶山羊背對著聲波收發設備時,拾音器在某一幀時刻下收到的回波中,包含呼吸和反芻運動的有效多徑信號較少,導致挑選不出符合規則的Si。

3.3 環境噪聲對監測結果的影響

與無干擾環境中的羊只相比,環境中存在其它羊叫聲干擾的Ⅵ號羊、音樂聲干擾的Ⅲ號羊和雨聲干擾的Ⅶ號羊,呼吸反芻監測漏幀率及呼吸反芻相對誤差并無明顯差別,如表4所示。

表4 靜臥奶山羊不同噪聲下的相對誤差和漏幀率

原因有兩方面。首先,與動物聲音特征分析的方法不同,本文是通過雷達原理,測量羊只呼吸和反芻過程中胸部及嘴部的周期性運動,完全不依賴于羊只自身發出的任何聲音。揚聲器發送特定聲波信號,聲波信號經過靜態環境及被測羊只發生一次或多次反射,最終被拾音器收到。另一方面,由于發送的超聲波頻段為38～40 kHz,遠超普通噪聲頻段,再結合通帶為38～40 kHz的帶通濾波器可以很好地在接收端屏蔽外界環境噪聲。因此系統對環境噪聲的抗干擾能力較強。

4 結論

(1)該方法可實現對單只靜臥奶山羊呼吸及反芻行為的同步監測。同時該方法對不同朝向的奶山羊魯棒,對環境噪聲具有較強的抗干擾能力。表明利用聲學沖激響應實現奶山羊呼吸及反芻的同步監測是可行的。

(2)該方法能準確地對靜臥奶山羊的呼吸及反芻行為進行監測。其中,呼吸平均相對誤差為2.60%,反芻平均相對誤差為3.51%,平均漏幀率為2.49%。