楊孫運，闞 秀

(上海工程技術大學 電子電氣工程學院,上海 201620)

近年來，隨著計算機圖像采集與處理技術的進步，各種監(jiān)測裝置的研發(fā)發(fā)展迅速。隨著視頻、圖像采集與監(jiān)測分辨率不斷提高，可利用現有檢測設備進行生物簡單活動的檢測跟蹤，并捕獲微小體態(tài)變化。在實驗室中，基于數字圖像處理獲取生物行為參數是目前生物行為研究中的熱點。行為參數是實驗生物行為分析的關鍵評價指標[1]。目前基于數字圖像處理的生物行為分析方法主要有基于閾值分割的生物行為分析方法、基于背景差分的生物行為分析方法、基于Meanshift跟蹤的生物行為分析方法和基于深度學習的生物行為分析方法[2-7]。

視覺傳感器拍攝到的圖像是二維的，只能檢測到生物在二維平面中運動的信息，無法獲得生物在三維空間的姿態(tài)信息。目前對于運動物體姿態(tài)信息的測量，多采用慣性測量傳感器(Inertial Measurement Unit，IMU)。九軸IMU是一款集成了三軸加速度計、三軸角速度計、三軸磁力計的微型裝置，具有體積小、重量輕、成本低、功耗低、可靠性高抗震動沖擊能力強等優(yōu)勢[8-10]，可以實時檢測出運動物體的姿態(tài)信息。鑒于此，在生物行為分析系統(tǒng)中，本文選用IMU測量小鼠的姿態(tài)信息。

本文利用數字圖像處理技術檢測出小鼠在二維平面上的運動行為信息；使用IMU傳感器檢測出小鼠在三維空間中的加速度、角速度、磁力等信息，并利用擴展卡爾曼濾波算法(Extended Kalman Filtering，EKF)解算出小鼠的姿態(tài)信息。實驗結果顯示，所設計的系統(tǒng)能夠有效地實現小鼠行為信息的全面檢測和分析。

1 系統(tǒng)綜述

圖1為所設計的系統(tǒng)的總體框圖，具體操作步驟如下：

圖1 視覺和IMU傳感器的生物行為分析系統(tǒng)圖Figure 1. Biological behavior analysis system diagram of vision and IMU sensors

步驟1將視覺傳感器固定在小鼠正上方進行小鼠活動拍攝，同時將IMU附著在小鼠的腹部進行參數測量；

步驟2對視頻中的圖像幀進行預處理。由于IMU姿態(tài)解算是實時的，所以想要得到小鼠的姿態(tài)信息，就需要設置響應指令記錄同一時刻IMU測量的參數；

步驟3利用主線程進行圖像幀的目標檢測和小鼠重心的計算，并利用子線程進行小鼠的姿態(tài)解算；

步驟4主線程和子線程全部結束后，獲得小鼠的重心和姿態(tài)角；

步驟5用當前幀的重心與上一幀的重心計算出位移、速度等運動行為參數，最后將運動行為參數和姿態(tài)參數保存并顯示到界面上。

循環(huán)上述步驟實現生物行為分析系統(tǒng)實時運行。

2 目標檢測及運動信息獲取

2.1 目標檢測

VIBE(Visual Background Extractor)算法是一種簡潔高效的運動目標檢測算法。該算法邏輯簡單，通過初始化背景模型、前景檢測、背景模型更新完成檢測任務[11]：

(1)背景模型初始化。VIBE算法利用視頻第一幀圖像初始化背景模型。在第一幀圖像中，任意一個像素值為v(x)的像素點x，在以其為中心的八鄰域中等概率隨機選取N個背景樣本元素作為像素點x的背景模型樣本集B(x)[12]，即

B(x)={B1(x),B2(x)，…,BN(x)}

(1)

式中，B(x)為像素點x的背景模型集合，以此完成背景模型初始化；

(2)前景檢測。在前景檢測過程中，將像素點x處新的像素值和背景樣本集中的元素比較判斷，確認是否屬于前景。樣本集合大小為N，以v(x)為中心，R為半徑建立一個二維歐式空間SR(v(x))，如圖2所示。若背景樣本集B(x)中的樣本落入SR(v(x))的采樣值個數大于給定的預設值Tmin，則認為v(x)是背景像素點，否則是前景像素點[13]；

圖2 二維歐式空間圖[14]Figure 2. Two-dimensional Euclidean space map

(3)背景模型更新。VIBE算法采用隨機更新策略，設定更新概率為φ。當前像素v(x)被判定為背景像素后，v(x)有1/φ的概率被更新為該像素點的背景樣本模型，隨機替換樣本模型中的任意一個樣本值，同時也有1/φ的概率隨機更新鄰域像素的背景模型的樣本值[14]。VIBE算法中，以上變量的默認值為N=20，R=30，Tmin=2，φ=16；

(4)算法更新。針對VIBE算法使用第一幀圖片初始化背景模型時，圖片中包含運動目標，在后續(xù)幀中易產生鬼影的問題。本文在實驗開始時通過加載無小鼠的初始背景幀解決鬼影問題，如圖3所示。

(a) (b) 圖3 改進前后的第20幀目標檢測圖(a)改進前 (b)改進后Figure 3.The 20th frame target detection map before and after the improvement(a)Before improvement (b)After improvement

針對實驗室內光線微弱變化的問題，本文使用自適應閾值抑制背景擾動。設η(x)為衡量光照變化影響的指標，則有

(2)

式中，ηi(x)表示第i幀的像素點x的像素值與背景模型樣本中像素值差分和的均值。ηi(x)越大，背景擾動越大；ηi(x)越小，背景越平穩(wěn)。匹配閾值Ri(x)的計算式如下

(3)

式中，α和δ為可控參數；Ri(x)表示像素點x在第i幀的匹配閾值。使用固定閾值和自適應閾值的結果對比如圖4。

(a) (b) 圖4 改進前后的第2 645幀目標檢測圖(a)改進前的第2 645幀目標檢測圖(b)改進后第2 645幀目標檢測圖Figure 4.The 2 645th frame target detection map before and after the improvement(a)The 2 645th frame target detection map before the improvement(b)The 2 645th frame target detection map after the improvement

2.2 運動行為信息計算

對檢測到的目標進行形態(tài)學處理，然后調用Opencv2中的findcontours函數獲得大鼠的輪廓點數組。隨后，調用moments函數對輪廓點數組處理以獲得大鼠輪廓的矩，并利用式(4)獲得輪廓的重心。

(4)

(5)

(6)

式中，Vn是速度，單位為：像素·s-1；tn是第n幀的時刻。

3 姿態(tài)解算

3.1 基礎知識

如圖5所示，導航坐標系為OXnYnZn，載體坐標系為OXbYbZb，其中Xn、Yn、Zn分別指向東、北、天，Xb、Yb、Zb分別指向小鼠的右、前、上。本文中，小鼠的姿態(tài)角是載體坐標系相對于導航坐標系的方向關系[15-16]。

圖5 導航坐標系和載體坐標系圖Figure 5. Navigation coordinate system and carrier coordinate system

導航坐標系下的參數轉換到載體坐標系下可以依次繞Z、Y、X旋轉φ、γ、θ得到。旋轉變換矩陣為式(7)。

(7)

由于四元數法可以避免歐拉角的奇異問題，計算量小且易于操作，也滿足實時性要求，因此系統(tǒng)采用四元數法來求解姿態(tài)角[16]。四元數是由4個元構成的數，如下式所示

q=q0+q1i+q2j+q3k

(8)

式中，q0、q1、q2、q3是實部，i、j、k為虛部。

導航坐標系到載體坐標系四元數形式的旋轉矩陣為

(9)

因此，將式(7)和式(9)聯(lián)立可解得四元數形式的姿態(tài)角。

(10)

3.2 基于EKF的姿態(tài)解算

3.2.1 狀態(tài)方程的建立

根據捷聯(lián)慣導系統(tǒng)的四元數理論[17]，四元數更新計算式如下

(11)

式中，Ts是采樣周期；k表示k時刻；wx、wy、wz為陀螺儀3個軸的測量角速度。故建立離散時間狀態(tài)方程為

Xk/k-1=Fk-1Xk-1/k-1+Wk-1

(12)

式中，Xk/k-1表示系統(tǒng)k-1時刻估計k時刻的狀態(tài)量，狀態(tài)量設為Xk=[q0(k)q1(k)q2(k)q3(k) ]T；Fk-1為k-1時刻的狀態(tài)轉移矩陣；Wk-1為k-1時刻系統(tǒng)的狀態(tài)噪聲矢量。狀態(tài)轉移矩陣Fk-1和過程噪聲方差矩陣Qk-1為如下形式

(13)

Qk-1=σqI4×4

(14)

式中，σq為姿態(tài)四元數系統(tǒng)噪聲；I4×4為4×4的單位矩陣。

3.2.2 觀測方程的建立

(15)

(16)

因此，本文將載體系下的加速度計和磁力計的測量值作為觀測量，觀測方程為

Zk=HkXk+Vk

(17)

式中，Zk表示k時刻的觀測量；Hk為k時刻的觀測矩陣；Vk為k時刻系統(tǒng)觀測噪聲。觀測量為Zk=

(18)

(19)

(20)

3.2.3基于EKF的姿態(tài)解算流程[18-19]

姿態(tài)解算具體步驟如下：

步驟1狀態(tài)量初始化。給定狀態(tài)量初值X0和誤差協(xié)方差矩陣P0；

步驟2狀態(tài)預測，如

Xk/k-1=Fk-1Xk-1/k-1

(21)

步驟3協(xié)方差預測，如

(22)

步驟4EKF增益，具體為

(23)

步驟5協(xié)方差更新，具體為

Pk/k=(I-KkHk)Pk/k-1

(24)

步驟6狀態(tài)更新，如下式

Xk/k=Xk/k-1+Kk[Zk-HkXk/k-1]

(25)

步驟7通過式(10)解算姿態(tài)角。

4 實驗與結果分析

4.1 實驗環(huán)境及實驗步驟

圖6 整體實驗環(huán)境圖(a)IMU傳感器 (b) 傳感器和小鼠 (c) 實驗平臺Figure 6. Overall experimental environment diagram(a)IMU sensor (b)Sensors and mice(c)The experiment platform

實驗步驟如下：

步驟1運行編寫好的軟件程序，出現生物行為分析系統(tǒng)界面，點擊控制區(qū)的加載按鍵加載實驗箱背景；

步驟2拖動控制區(qū)的行列滑動塊設置周邊區(qū)域和中間區(qū)域的大?。?/p>

步驟3將附著好傳感器的小鼠放入實驗箱，點擊開始按鍵；

步驟4系統(tǒng)運行9 min后，點擊退出按鍵。

4.2 結果分析

如圖7所示，實驗開始后，小鼠在實驗平臺中的自發(fā)活動可以通過生物行為分析系統(tǒng)的界面實時顯示。控制區(qū)可以控制實驗平臺中的周邊區(qū)域和中間區(qū)域的大小，例如原視頻中的矩形框。參數顯示區(qū)顯示小鼠當前時刻的行為參數信息。

圖7 實時結果顯示圖Figure 7. Real-time results display diagram

系統(tǒng)運行結束后，輸出小鼠在整個實驗中的行為分析結果。統(tǒng)計后如表1和圖8～圖10所示。

表1 小鼠實驗結果參數表Table 1. Parameters of experimental results in mice

圖8 俯仰角Figure 8. Pitching angle

圖9 滾轉角Figure 9. Roll angle

圖10 偏航角Figure 10. Yaw angle

由表1可以看出小鼠在進入實驗箱初期，活動量較少，大部分時間停留在周邊或墻角；隨著時間增加，小鼠對于試驗箱環(huán)境逐漸熟悉，活動量增加，在實驗箱中央區(qū)域的時間也逐漸增加，實驗結果與所有曠場實驗結果一致。

由圖8可以看出小鼠在9分鐘實驗中的俯仰角變化情況，當俯仰角先增大后減小且極大值大于30°時認為發(fā)生了站立姿態(tài)。在誤差允許的范圍內，此次實驗小鼠站立次數為8次。綜合圖8和圖9，可以看出小鼠在站立以外的其他時刻俯仰角和滾轉角在-20°～ +20°之間波動。圖10反映了小鼠在試驗箱中的方向變化情況。

5 結束語

實驗室內小鼠行為分析是生物行為分析領域重要的研究課題之一。本文在Python編程環(huán)境下，提出利用數字圖像處理技術并結合IMU傳感器檢測小鼠的運動行為信息和姿態(tài)信息。實驗結果表明，該系統(tǒng)可以較好地運行，并實時檢測出小鼠的運動行為信息和姿態(tài)信息。本文所設計的生物行為分析系統(tǒng)相較于傳統(tǒng)人工觀察法，可節(jié)約成本，數據可靠。相較于單一的圖像處理方法，本文方法既獲得了實驗生物的運動行為信息，又獲得了實驗生物的姿態(tài)信息。這些行為信息可被用來研究實驗動物的神經功能、心里過程、藥物作用等，為神經科學、認知科學、藥物篩選以及行為和基因的相關性研究提供了參考。