基于图像和IMU传感器的生物行为分析系统设计

杨孙运，阚 秀

(上海工程技术大学 电子电气工程学院,上海 201620)

近年来，随着计算机图像采集与处理技术的进步，各种监测装置的研发发展迅速。随着视频、图像采集与监测分辨率不断提高，可利用现有检测设备进行生物简单活动的检测跟踪，并捕获微小体态变化。在实验室中，基于数字图像处理获取生物行为参数是目前生物行为研究中的热点。行为参数是实验生物行为分析的关键评价指标[1]。目前基于数字图像处理的生物行为分析方法主要有基于阈值分割的生物行为分析方法、基于背景差分的生物行为分析方法、基于Meanshift跟踪的生物行为分析方法和基于深度学习的生物行为分析方法[2-7]。

视觉传感器拍摄到的图像是二维的，只能检测到生物在二维平面中运动的信息，无法获得生物在三维空间的姿态信息。目前对于运动物体姿态信息的测量，多采用惯性测量传感器(Inertial Measurement Unit，IMU)。九轴IMU是一款集成了三轴加速度计、三轴角速度计、三轴磁力计的微型装置，具有体积小、重量轻、成本低、功耗低、可靠性高抗震动冲击能力强等优势[8-10]，可以实时检测出运动物体的姿态信息。鉴于此，在生物行为分析系统中，本文选用IMU测量小鼠的姿态信息。

本文利用数字图像处理技术检测出小鼠在二维平面上的运动行为信息；使用IMU传感器检测出小鼠在三维空间中的加速度、角速度、磁力等信息，并利用扩展卡尔曼滤波算法(Extended Kalman Filtering，EKF)解算出小鼠的姿态信息。实验结果显示，所设计的系统能够有效地实现小鼠行为信息的全面检测和分析。

1 系统综述

图1为所设计的系统的总体框图，具体操作步骤如下：

图1 视觉和IMU传感器的生物行为分析系统图Figure 1. Biological behavior analysis system diagram of vision and IMU sensors

步骤1将视觉传感器固定在小鼠正上方进行小鼠活动拍摄，同时将IMU附着在小鼠的腹部进行参数测量；

步骤2对视频中的图像帧进行预处理。由于IMU姿态解算是实时的，所以想要得到小鼠的姿态信息，就需要设置响应指令记录同一时刻IMU测量的参数；

步骤3利用主线程进行图像帧的目标检测和小鼠重心的计算，并利用子线程进行小鼠的姿态解算；

步骤4主线程和子线程全部结束后，获得小鼠的重心和姿态角；

步骤5用当前帧的重心与上一帧的重心计算出位移、速度等运动行为参数，最后将运动行为参数和姿态参数保存并显示到界面上。

循环上述步骤实现生物行为分析系统实时运行。

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算法使用第一帧图片初始化背景模型时，图片中包含运动目标，在后续帧中易产生鬼影的问题。本文在实验开始时通过加载无小鼠的初始背景帧解决鬼影问题，如图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)越小，背景越平稳。匹配阈值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 运动行为信息计算

对检测到的目标进行形态学处理，然后调用Opencv2中的findcontours函数获得大鼠的轮廓点数组。随后，调用moments函数对轮廓点数组处理以获得大鼠轮廓的矩，并利用式(4)获得轮廓的重心。

(4)

(5)

(6)

式中，Vn是速度，单位为：像素·s-1；tn是第n帧的时刻。

3 姿态解算

3.1 基础知识

如图5所示，导航坐标系为OXnYnZn，载体坐标系为OXbYbZb，其中Xn、Yn、Zn分别指向东、北、天，Xb、Yb、Zb分别指向小鼠的右、前、上。本文中，小鼠的姿态角是载体坐标系相对于导航坐标系的方向关系[15-16]。

图5 导航坐标系和载体坐标系图Figure 5. Navigation coordinate system and carrier coordinate system

导航坐标系下的参数转换到载体坐标系下可以依次绕Z、Y、X旋转φ、γ、θ得到。旋转变换矩阵为式(7)。

(7)

由于四元数法可以避免欧拉角的奇异问题，计算量小且易于操作，也满足实时性要求，因此系统采用四元数法来求解姿态角[16]。四元数是由4个元构成的数，如下式所示

q=q0+q1i+q2j+q3k

(8)

式中，q0、q1、q2、q3是实部，i、j、k为虚部。

导航坐标系到载体坐标系四元数形式的旋转矩阵为

(9)

因此，将式(7)和式(9)联立可解得四元数形式的姿态角。

(10)

3.2 基于EKF的姿态解算

3.2.1 状态方程的建立

根据捷联惯导系统的四元数理论[17]，四元数更新计算式如下

(11)

式中，Ts是采样周期；k表示k时刻；wx、wy、wz为陀螺仪3个轴的测量角速度。故建立离散时间状态方程为

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

(12)

式中，Xk/k-1表示系统k-1时刻估计k时刻的状态量，状态量设为Xk=[q0(k)q1(k)q2(k)q3(k) ]T；Fk-1为k-1时刻的状态转移矩阵；Wk-1为k-1时刻系统的状态噪声矢量。状态转移矩阵Fk-1和过程噪声方差矩阵Qk-1为如下形式

(13)

Qk-1=σqI4×4

(14)

式中，σq为姿态四元数系统噪声；I4×4为4×4的单位矩阵。

3.2.2 观测方程的建立

(15)

(16)

因此，本文将载体系下的加速度计和磁力计的测量值作为观测量，观测方程为

Zk=HkXk+Vk

(17)

式中，Zk表示k时刻的观测量；Hk为k时刻的观测矩阵；Vk为k时刻系统观测噪声。观测量为Zk=

(18)

(19)

(20)

3.2.3基于EKF的姿态解算流程[18-19]

姿态解算具体步骤如下：

步骤1状态量初始化。给定状态量初值X0和误差协方差矩阵P0；

步骤2状态预测，如

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

(21)

步骤3协方差预测，如

(22)

步骤4EKF增益，具体为

(23)

步骤5协方差更新，具体为

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

(24)

步骤6状态更新，如下式

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

(25)

步骤7通过式(10)解算姿态角。

4 实验与结果分析

4.1 实验环境及实验步骤

图6 整体实验环境图(a)IMU传感器 (b) 传感器和小鼠 (c) 实验平台Figure 6. Overall experimental environment diagram(a)IMU sensor (b)Sensors and mice(c)The experiment platform

实验步骤如下：

步骤1运行编写好的软件程序，出现生物行为分析系统界面，点击控制区的加载按键加载实验箱背景；

步骤2拖动控制区的行列滑动块设置周边区域和中间区域的大小；

步骤3将附着好传感器的小鼠放入实验箱，点击开始按键；

步骤4系统运行9 min后，点击退出按键。

4.2 结果分析

如图7所示，实验开始后，小鼠在实验平台中的自发活动可以通过生物行为分析系统的界面实时显示。控制区可以控制实验平台中的周边区域和中间区域的大小，例如原视频中的矩形框。参数显示区显示小鼠当前时刻的行为参数信息。

图7 实时结果显示图Figure 7. Real-time results display diagram

系统运行结束后，输出小鼠在整个实验中的行为分析结果。统计后如表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可以看出小鼠在进入实验箱初期，活动量较少，大部分时间停留在周边或墙角；随着时间增加，小鼠对于试验箱环境逐渐熟悉，活动量增加，在实验箱中央区域的时间也逐渐增加，实验结果与所有旷场实验结果一致。

由图8可以看出小鼠在9分钟实验中的俯仰角变化情况，当俯仰角先增大后减小且极大值大于30°时认为发生了站立姿态。在误差允许的范围内，此次实验小鼠站立次数为8次。综合图8和图9，可以看出小鼠在站立以外的其他时刻俯仰角和滚转角在-20°～ +20°之间波动。图10反映了小鼠在试验箱中的方向变化情况。

5 结束语

实验室内小鼠行为分析是生物行为分析领域重要的研究课题之一。本文在Python编程环境下，提出利用数字图像处理技术并结合IMU传感器检测小鼠的运动行为信息和姿态信息。实验结果表明，该系统可以较好地运行，并实时检测出小鼠的运动行为信息和姿态信息。本文所设计的生物行为分析系统相较于传统人工观察法，可节约成本，数据可靠。相较于单一的图像处理方法，本文方法既获得了实验生物的运动行为信息，又获得了实验生物的姿态信息。这些行为信息可被用来研究实验动物的神经功能、心里过程、药物作用等，为神经科学、认知科学、药物筛选以及行为和基因的相关性研究提供了参考。