CN108765484B - 基于两台高速摄像机的活体昆虫运动采集及数据重构方法 - Google Patents

Abstract

本公开提供了一种基于两台高速摄像机的活体昆虫运动采集及数据重构方法，两台高速摄像机的光轴垂直于供活体昆虫活动的观察箱的壁面，并且两台高速摄像机的光轴正交；该方法包括：建立高速摄像机的坐标变换模型；通过基于共面点的分步线性标定方法求解高速摄像机的参数；分别获取观测点在两台高速摄像机中的场深；以及在两台高速摄像机所拍摄的图像中分别提取具有相同特征点的图像坐标来完成对应点的匹配，并根据高速摄像机参数和场深对观测点进行三维重构。

Description

基于两台高速摄像机的活体昆虫运动采集及数据重构方法

技术领域

本公开涉及一种基于两台高速摄像机的活体昆虫运动采集及数据重构方法。

背景技术

针对昆虫飞行的空气动力学研究和飞行动力学研究都需要获得昆虫扑翼飞行的真实运动规律，因此，通过高速摄像机记录飞行昆虫影像，处理分析得到昆虫扑翼飞行中身体和翅膀的运动学参数十分必要。

起初由于设备限制，研究者多利用食物、气味或是利用其固有的趋光性引导其飞行。1997年，Willmott和Ellington利用喂食器令鹰蛾悬停在拍摄画面内，从而完成了对鹰蛾的扑翼频率和翅膀扭转角的测量。虽然通过引导昆虫飞行可以获得接近真实飞行的运动规律，但只能构造出悬停飞行或直线飞行等有限的飞行状态，无法满足空气动力学以及飞行动力学研究中对于运动学参数多样性的要求。因此，发展能够实现对昆虫自由飞行进行运动学参数测量的实验方法势在必行。

通常的测量方法是，先确定一个适当大小的拍摄区域，当昆虫进入该区域后即触发摄像机进行拍摄，从而捕获昆虫经过该区域的飞行图像。这种测量方法的问题是，若拍摄区域过大，则昆虫所占像素较少，分辨率无法达到测量要求；若拍摄区域过小，则昆虫进入拍摄区域的几率以及在区域内停留的时间也相应减少，降低了实验的成功率。而且，对于自由飞行的观测实验，昆虫在拍摄区域内的身体姿态和运动轨迹都是不确定的，所以运动学参数的测量需要通过后期的三维重构完成。

Azuma在1988年通过控制风洞的风速使蜻蜓呈现出近似悬停的低速前飞状态，利用实验前在蜻蜓翅膀上绘制的线型标记，测量获得了蜻蜓飞行状态下的翅膀攻角和拍动角。这一实验仅借助一台胶片式高速摄像机完成。Sato和Azuma又采用类似的方法测量了豆娘飞行的运动学参数。王浩、曾理江等人利用梳状条纹投影成功测量了蜻蜓自由飞行中翅膀的运动学参数以及翅膀的变形。在摄像机拍摄区域内投射梳状激光条纹，通过高速摄像机记录条纹在蜻蜓翅膀上投影的位置和形状，之后通过插值求解其运动和变形的相关参数。上述方法主要是利用1台摄像机拍摄，结合相应的光学方法及假设进行三维数据重构，但是这种方法存在两点不足：一是对于实验状态要求苛刻，实验成功率较低；二是为了获得更高的分辨率，拍摄区域均较小，因而无法获得长时间连续的飞行图像。

1997年Wakeling和Ellington在剑桥大学内搭建了一个大型温室，温室内包括一个8.9mx1.3mx0.7m的池塘，池塘周围按照真实的蜻蜓栖息地的生态环境种植了相应的水草和灌木丛。通过两台光轴正交的高速摄像机拍摄蜻蜓自由飞行。三维重构采用三角测量法，在摄像机视野内悬挂有三根不同长度的标记棒，根据图像中蜻蜓与标记棒的相对位置关系即可求得蜻蜓的身体和翅膀的姿态，从而获得其运动学参数。通过这一方法，Wakeling和Ellington测量了蜻蜓和豆娘自由飞行时的运动学参数，并对两种昆虫的飞行进行了对比。该方法的问题是试验成本高，不利于推广，同时三维重构方法并不能实现翅膀精确角度(攻角、偏离角等)参数的测量。

Fry和Sun通过三台光轴彼此正交的高速红外摄像机测量了果蝇和食蚜蝇机动转弯飞行时身体位移、速度和翅膀拍动角、攻角等运动学参数。在后期三维重构过程中，根据真实昆虫的几何信息构建了身体和翅膀的简化模型，并且在视场中不断调整其位置和姿态，直到模型在三个视角的投影均与相应图像吻合。该方法是基于三台摄像机，试验成本较高；三维重构方法上存在两点不足：一是搭建的数学模型需要不断调整以便图像拍摄，时间成本比较高，同时精度并不能很好的保证；二是为了获得更高的分辨率，采用的拍摄区域均较小，因而无法获得长时间连续的飞行图像。

综上所述，目前活体昆虫拍摄及数据重构技术主要存在以下问题：(1)长时间连续的飞行图像的获取与三维数据重构方法的简易性和精确性；(2)实验方法的简易性和成本的控制。

摄像机标定是三维数据重构必不可少的一个步骤，关于摄像机的标定方法，邹凤娇提出了一种简单快速、标定精度高的基于共面点的分步线性标定方法。该方法通过分步标定并建立一种新的畸变模型,实现了只利用共面点标定物,不需要摄像机做任何运动,就可以全部线性地求解出摄像机的内外参数。既避免了非线性优化的繁琐和不稳定,也解决了其它线性法中部分内参数的标定问题，值得借鉴参考。

发明内容

为了解决至少一个上述技术问题，本公开提供了一种基于两台高速摄像机的活体昆虫运动采集及数据重构方法。

根据本公开的一个方面，基于两台高速摄像机的活体昆虫运动采集及数据重构方法中，

两台高速摄像机的光轴垂直于供活体昆虫活动的观察箱的壁面，并且两台高速摄像机的光轴正交；

该方法包括：

建立高速摄像机的坐标变换模型，其中，坐标变换模型包括：刚体变换模型，用于从世界坐标系到摄像机坐标系的坐标变换；透视投影模型，用于从摄像机坐标系到理想成像平面坐标系的坐标变换；径向畸变模型，用于从理想成像平面坐标系到实际成像平面坐标系的坐标变换；以及像素单位变换，用于从实际成像平面坐标系到图像坐标系的坐标变换；

通过建立坐标变换模型来确定需要求解的高速摄像机的参数；

通过基于共面点的分步线性标定方法求解高速摄像机的参数；

分别获取该观测点在两台高速摄像机中的场深；以及

在两台高速摄像机所拍摄的图像中分别提取具有相同特征点的图像坐标来完成对应点的匹配，并根据高速摄像机的参数和场深对该观测点进行三维重构。

根据本公开的至少一个实施方式，

刚体变换模型由一个旋转矩阵R和一个平移向量T来描述，相应齐次坐标之间的关系如下式1和2：

其中，(x_c,y_c,z_c)表示观测点在高速摄像机的坐标系下的三维坐标，(X_w,Y_w,Z_w)表示观测点在世界坐标系下的三维坐标；

透视投影模型的相应齐次坐标之间的关系如下式3：

其中，f为高速摄像机的有效焦距，(x_u,y_u)表示观测点在成像平面坐标系下的理想成像坐标；

径向畸变模型是在透视投影模型的基础上加入畸变的影响，径向畸变模型仅考虑三阶径向畸变，相应齐次坐标之间的关系如下式4：

其中，(x_d,y_d)表示观测点在成像平面坐标系下的实际成像坐标，k_i表示各阶径向畸变系数，r表示像点到高速摄像机的光轴的距离；

像素单位变换是由物理长度向像素单位的变换，实际成像平面坐标与图像坐标之间的关系如下式5：

其中，(u,v)表示实际成像坐标在图像坐标系下的坐标，(u₀,v₀)表示成像平面坐标系原点在图像坐标系下的坐标，dx表示单位图像元素在水平方向的尺寸大小，dy表示单位图像元素在竖直方向的尺寸大小，

表示图像元素比例因子。

根据本公开的至少一个实施方式，摄像机参数包括：

摄像机内部参数，包括成像平面坐标系的原点在图像坐标系下的坐标(u₀,v0)、图像元素比例因子

有效焦距f和径向畸变系数k_i；以及

摄像机外部参数，包括旋转矩阵R和平移向量T。

根据本公开的至少一个实施方式，摄像机参数的求解方法包括：

将标定所用的世界坐标系OX_wY_wZ_w固定在二维标定板上，世界坐标系的X_w轴和Y_w轴平行于标定板，坐标系符合右手法则；

将二维标定板放置在高速摄像机的拍摄空间中，提取二维标定板的图像中各标定点的图像坐标(u,v)；

根据式5求解观测点在成像平面坐标系下的实际成像坐标(x_d,y_d)；

根据径向约束条件：

获得下式6：

将式6右端项的分子和分母同时乘以1/t_y，即将透镜变形对有效焦距f的影响考虑进来，获得下式7：

采用最小二乘方法求解式7的超定方程组并获得5个中间变量，中间变量分别为：

根据I₁、I₂、I₄、I₅求解t_y，如果I₁、I₂、I₄、I₅之中的任意两项为零，则

等于两个非零项平方和的倒数，否则

的值由下式9计算：

其中，

B＝(I₁I₅-I₂I₄)²；

假设t_y为正，由式8分别求出r₁₁、r₁₂、t_x、r₂₁和r₂₂的值，并将r₁₁、r₁₂、t_x、r₂₁和r₂₂的值代入式6计算得到x_c和y_c；

根据旋转矩阵R的正交性继续求解R的其它元素，首先根据下式10计算r₁₃和r₂₃：

此处同样先假设r₁₃符号为正，而r₂₃的符号则与r₁₁r₂₁+r₁₂r₂₂的符号相反，旋转矩阵R第三行的元素分别为：

联立式3和式4可以得到下式12：

将上式12整理为矩阵形式，如下式13：

其中，E＝r₁₁X_w+r₁₂Y_w+t_x，F＝r₂₁X_w+r₂₂Y_w+t_y，G＝r₃₁X_w+r₃₂Y_w，采用最小二乘方法求解上式13的超定方程组即可获得有效焦距f、径向畸变系数k_i和平移分量t_z。

根据本公开的至少一个实施方式，

假设t_y为正，由式8分别求出r₁₁、r₁₂、t_x、r₂₁和r₂₂的值，并将r₁₁、r₁₂、t_x、r₂₁和r₂₂的值代入式6计算得到x_c和y_c；如果x_c和y_c的符号与x_d和y_d的符号一致，则证明假设是正确的；

如果x_c和y_c的符号与x_d和y_d不一致，则将t_y的符号改为负，重新计算r₁₁、r₁₂、t_x、r₂₁和r₂₂的值，并将新的r₁₁、r₁₂、t_x、r₂₁和r₂₂的值代入式6重新计算x_c和y_c。

根据本公开的至少一个实施方式，如果求解得到的有效焦距f为负数，则回到式10将r₁₃和r₂₃改为相反的符号，并根据式11重新计算旋转矩阵R的第三行元素，进而更新式13的系数矩阵，即可获得正确的有效焦距f、径向畸变系数k_i和平移分量t_z。

根据本公开的至少一个实施方式，获取场深的方法包括：

将两台高速摄像机呈光轴正交放置；以及

观测点在一台高速摄像机的图像中沿水平方向的距离是观测点在另一台摄像机中的场深。

根据本公开的至少一个实施方式，对高速摄像机进行标定之前，利用圆环标定板对两台高速摄像机的光轴的垂直度进行校准。

根据本公开的至少一个实施方式，

对高速摄像机进行标定所借助的标定物为二维棋盘标定板；

二维棋盘标定板的尺寸为600mmx900mm，棋盘图案中每一格的尺寸为100mmx100mm；以及

二维棋盘标定板相对于高速摄像机的光轴倾斜放置在观察箱中。

根据本公开的至少一个实施方式，高速摄像机的拍摄区域设置在观察箱的中间区域。

附图说明

附图示出了本公开的示例性实施方式，并与其说明一起用于解释本公开的原理，其中包括了这些附图以提供对本公开的进一步理解，并且附图包括在本说明书中并构成本说明书的一部分。

图1是根据本公开的至少一个实施方式的蜻蜓扑翼飞行实验观测系统的组成及布局示意图。

图2是根据本公开的至少一个实施方式的摄像机成像过程中四个坐标系的示意图。

图3是根据本公开的至少一个实施方式的坐标系之间的坐标变换关系示意图。

图4是根据本公开的至少一个实施方式的透视投影模型原理示意图。

图5是根据本公开的至少一个实施方式的观察箱壁面上的标尺示意图。

图6是根据本公开的至少一个实施方式的观察箱上的世界坐标系示意图。

图7是根据本公开的至少一个实施方式的摄像机光轴与观察箱壁面存在夹角时的成像示意图。

图8是根据本公开的至少一个实施方式的用于光轴垂直度校准的圆环标定板。

图9是根据本公开的至少一个实施方式的Image Pro Plus6.0工作界面示意图。

图10是根据本公开的至少一个实施方式的摄像机标定所用二维标定板示意图。

图11是根据本公开的至少一个实施方式的两台摄像机拍摄的标定图像。

图12是根据本公开的至少一个实施方式的二维标定板的拍摄图像与坐标提取示意图。

图13是根据本公开的至少一个实施方式的二维标定板上20个标定点的三维重构结果示意图。

图14是根据本公开的至少一个实施方式的蜻蜓起飞过程翅膀运动参数图。

图15是根据本公开的至少一个实施方式的蜻蜓起飞过程质心运动图。

具体实施方式

下面结合附图和实施方式对本公开作进一步的详细说明。可以理解的是，此处所描述的具体实施方式仅用于解释相关内容，而非对本公开的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本公开相关的部分。

需要说明的是，在不冲突的情况下，本公开中的实施方式及实施方式中的特征可以相互组合。下面将参考附图并结合实施方式来详细说明本公开。

开展昆虫扑翼飞行观测实验的最终目的是，基于计算机视觉技术，根据摄像机所拍摄的图像获取昆虫在三维空间内的飞行轨迹与飞行姿态。为此，则必须建立三维空间与二维图像之间的对应关系。

实验分析所用图像由摄像机拍摄所得，其与三维空间的对应关系完全由摄像机成像模型以及摄像机在空间中的位置决定。考虑到昆虫扑翼飞行观测实验精度要求较高，而且完成标定后摄像机的方位和相关设置在拍摄过程中保持不变。因此，在本公开的一个可选实施方式中，选用需要借助标定物的传统标定方法对实验所用的两台高速摄像机进行标定。所标定参数的具体求解方法选用基于共面点的分步线性标定方法。该方法利用共面点约束，只需借助二维标定物即可完成标定，操作简单，容易实现。通过分步思想，避免了求解的非线性迭代，仅通过线性解析计算即可获得全部参数，节省了大量计算资源，有效提高了标定速度。同时，在标定中考虑了必要的径向畸变以及透镜变形等因素的影响，故而得以保证较高的精度。

下面结合附图并通过蜻蜓扑翼飞行实验来详细说明活体昆虫的运动采集及数据重构方法。

在本公开的一个可选实施方式中，如图1所示，为蜻蜓扑翼飞行实验的观测系统，包括两台高速摄像机、与高速摄像机配套的控制器显示装置(CDU)、两盏摄影灯和一个观察箱。

高速摄像机的选择主要取决于其录制速度(即帧速率)以及在相应录制速度下所能保证的最大图像分辨率。本公开选用的两台高速摄像机均为Olympus i-SPEED TR高速摄像机。蜻蜓扑翼飞行过程中翅膀的运动速度较高，很有可能在一帧周期内发生可察觉的位移，为了获取清晰地飞行图像，避免拖影等现象的产生，通常应该适当缩短快门时间，并采用两盏功率为1000W的摄影灯为拍摄区域提供照明补偿。观察箱的尺寸为780mmx780mmx1000mm，观察箱壁面为玻璃材质。两台高速摄像机的光轴呈正交状态，可以转化为两台高速摄像机的光轴分别垂直于正对的观察箱的一侧壁面。在本公开的一个可选实施方式中，如图1所示，当蜻蜓进入观察箱后，一般会主动开始飞行，大部分蜻蜓会直接飞向观察箱顶部，由于靠近观察箱顶部和底部的飞行并不能很好的表征蜻蜓自由飞行的特性，因此，还需要将两台摄像机的拍摄区域选定在观察箱的中间区域。

首先，建立高速摄像机的坐标变换模型，其中，该坐标变换模型包括：刚体变换模型，用于从世界坐标系到摄像机坐标系的坐标变换；透视投影模型，用于从摄像机坐标系到理想成像平面坐标系的坐标变换；径向畸变模型，用于从理想成像平面坐标系到实际成像平面坐标系的坐标变换；以及像素单位变换，用于从实际成像平面坐标系到图像坐标系的坐标变换。

在完成上述坐标变换模型的建立之后，高速摄像机标定过程中需要求解的全部内外参数也随之确定。通过基于共面点的分步线性标定方法求解高速摄像机的全部内外参数，即可以获取活体蜻蜓身体上某一观测点的三维坐标与所述观测点对应的图像坐标之间的关系。再分别获取观测点在两台高速摄像机中的场深，通过从两台高速摄像机所拍摄的图像中提取相同特征点的图像坐标来完成对应点的匹配，从而对观测点进行三维重构。

在本公开的一个可选实施方式中，高速摄像机的完整成像过程如图2所示，从观测点的世界坐标P(X_w,Y_w,Z_w)到观测点对应的图像坐标p(u,v)，共需要经过四次坐标变换。图中的四个坐标系分别为世界坐标系OX_wY_wZ_w，摄像机坐标系Ox_cy_cz_c，成像平面坐标系O₁xy和图像坐标系O₀uv。

世界坐标系OX_wY_wZ_w：世界坐标系为三维坐标系，可以根据实际情况任意选取，是一个人为确定的坐标系，观测点P在世界坐标系下的三维坐标可以记为(X_w,Y_w,Z_w)。

摄像机坐标系Ox_cy_cz_c：摄像机坐标系为三维坐标系，是固定在高速摄像机上的坐标系。摄像机坐标系的原点O位于摄像机的光心，z_c轴与摄像机光轴重合，方向指向摄像机前方，x_c轴和y_c轴平行于摄像机成像平面，其方向满足右手法则。观测点P在摄像机坐标系下的三维坐标可以记为(x_c,y_c,z_c)。

成像平面坐标系O₁xy：成像平面坐标系为二维坐标系，是固定在摄像机成像平面上的坐标系。成像平面坐标系的原点O₁位于摄像机光轴与成像平面的交点处，x轴和y轴均位于摄像机成像平面内，x轴和y轴分别平行于x_c和y_c轴。摄像机光心到成像平面的距离即

为摄像机的有效焦距f。观测点P经过摄像机后的理想成像位置P_u和实际成像位置P_d均位于成像平面内，二者在成像平面坐标系下的二维坐标分别记为P_u(x_u,y_u)与P_d(x_d,y_d)。

图像坐标系O₀uv：图像坐标系为二维坐标系，是固定在数字图像上的坐标系。与上述三个坐标系不同的是，图像坐标系并不是以物理长度为单位的，而是以图像元素为单位。图像坐标系的原点O₀位于数字图像的左上角，u轴和v轴分别对应图像的横向和纵向，u轴正方向向右，v轴正方向向下。观测点P的实际成像位置P_d在图像坐标系下的二维坐标可以记为(u,v)。该坐标表示的分别是这一像素点在图像中的列数和行数。成像平面坐标系的原点O₁在图像坐标系下的坐标可以记为(u₀,v₀)，O₁被称为图像主点，(u₀,v₀)称为图像主点坐标，理论上图像主点位于图像的中心。

从世界坐标系到图像坐标系，需先后经历刚体变换模型、透视投影模型、径向畸变模型以及像素单位变换四次坐标变换，综合以上四次变换的模型称为摄像机的坐标变换模型。上述各坐标系之间的坐标变换关系如图3所示。下面结合附图2、图3和图4详细说明高速摄像机的坐标变换模型。

(1)刚体变换模型:从世界坐标系到摄像机坐标系的坐标变换是刚体变换，可以由一个旋转矩阵R和一个平移向量T描述，相应齐次坐标之间的关系如下式1和2：

其中，(x_c,y_c,z_c)表示观测点在高速摄像机的坐标系下的三维坐标，(X_w,Y_w,Z_w)表示观测点在世界坐标系下的三维坐标。

(2)透视投影模型：从摄像机坐标系到理想成像平面坐标系的坐标变换满足透视投影模型。透视投影模型即摄像机的理想模型。如图4所示，在三维空间中的物体上的一点(物点)记为Q，Q点在摄像机坐标系Oxyz下的坐标为(X,Y,Z)，Q点在成像平面上的像点q在摄像机坐标系Oxyz下的坐标为(x,y,f)。其中f为成像模型的焦距。由针孔成像原理可知，物点Q在成像平面上的成像位置q即位于透镜中心O与物点Q的连线OQ与成像平面的交点。因此，在摄像机的理想模型中，像点q与物点Q的坐标之间的关系可表示为下式14：

由式14可获得透视投影模型的相应齐次坐标之间的关系，如下式3：

其中，f为高速摄像机的有效焦距，(x_u,y_u)表示观测点在成像平面坐标系下的理想成像坐标。

(3)径向畸变模型：从理想成像平面坐标系到实际成像平面坐标系的坐标变换，即在摄像机理想模型的基础上加入畸变的影响。

由于摄像机成像系统中的光学组件在制作和装配过程中必然存在误差，所以实际成像过程并不是理想的透视投影模型。若要获得较高的标定精度，需要在摄像机模型建立时考虑实际成像的各类畸变。根据畸变产生的原因不同，实际成像中存在的光学畸变可分为径向畸变、偏心畸变和薄棱镜畸变三大类。但是建立摄像机模型时并不需要考虑全部的非线性因素，过于复杂的模型将会在摄像机标定过程中引入大量非线性参数，不仅不能达到提高精度的目的，反而会适得其反，造成算法的不稳定。一般来讲，对于普通镜头只需要针对径向畸变予以修正即可满足测量要求。因此，在本公开中高速摄像机标定过程采用仅考虑三阶径向畸变的畸变模型。

径向畸变是由于透镜实际的几何形状与理论形状之间存在误差而产生的，属于对称畸变。由于径向畸变的存在，实际成像位置相对于理想成像位置会产生径向偏移。如果像点沿径向靠近光轴，则径向畸变为负，称为桶形畸变；反之，如果像点沿径向远离光轴，则径向畸变为正，称为枕形畸变。径向畸变是摄像机实际成像畸变中的主要分量，径向畸变的一般形式可表示为下式15：

其中，k_i表示各阶径向畸变系数，

表示像点到高速摄像机的光轴的距离。可见越靠近图像边缘径向畸变越严重。

考虑三阶径向畸变的畸变模型，观测点的实际成像位置(x_d,y_d)与理想成像位置(x_u,y_u)之间存在如下式16的关系：

由式16可获得径向畸变模型的相应齐次坐标之间的关系，如下式4：

其中，(x_d,y_d)表示观测点在成像平面坐标系下的实际成像坐标。

像素单位变换：从实际成像平面坐标系到图像坐标系的坐标变换是由物理长度到像素单位的变换。观测点的实际成像坐标(x_d,y_d)与图像坐标(u，v)之间的关系如下式5：

其中，(u₀,v₀)表示成像平面坐标系原点(图像主点O₁)在图像坐标系下的坐标，dx表示单位图像元素在水平方向的尺寸大小，dy表示单位图像元素在竖直方向的尺寸大小，

表示图像元素比例因子。

在本公开的一个可选实施方式中，完成上述高速摄像机坐标变换模型的建立之后，高速摄像机标定过程中需要求解的全部内外参数也随之确定。包括图像主点坐标(u₀,v₀)、图像元素比例因子s_x、有效焦距f和径向畸变系数k_i等由摄像机内部光学组件性质决定的内部参数，以及旋转矩阵R和平移向量T等表征摄像机在世界坐标系中的具体位置及方位的外部参数。

在本公开的一个可选实施方式中，高速摄像机的标定过程即求解上述内外参数的过程。本公开采用的求解方法是基于共面点的分步线性标定方法，该方法的关键在于应用了径向约束条件，该径向约束条件的数学表达式如下式17：

将标定所用的世界坐标系OX_wY_wZ_w固定在二维标定板上。世界坐标系的X_w轴和Y_w轴平行于二维标定板，即所有标定点的Z_w坐标为零，坐标系符合右手法则。

图像主点O₁位于图像中心，本实施方式所用高速摄像机的分辨率为1280x1024。因此图像主点O₁(u₀,v₀)的坐标分别为u₀＝640，v₀＝512。根据高速摄像机说明书中的标注，其单位图像元素尺寸dx和dy均为0.021mm，则图像元素比例因子s_x＝1。

首先，通过高速摄像机拍摄在观察箱中任意放置的二维标定板，并提取图像中各标定点的图像坐标(u,v)。

然后，根据已知的上述高速摄像机的内部参数，以及上式5，求解观测点在成像平面坐标系下的实际成像坐标(x_d,y_d)。

由式17可以获得下式6：

将式6右端项的分子和分母同时乘以1/t_y，即将透镜变形对有效焦距f的影响考虑进来，稍加整理后，可以获得下式7：

上式7是关于5个中间变量的线性方程，一般来讲，落在高速摄像机视野内的标定点的个数都远远大于5个，因此采用最小二乘方法求解式7的超定方程组，并获得5个中间变量，中间变量分别为：

因为旋转矩阵R是单位正交矩阵，因此根据I₁、I₂、I₄、I₅之间的正交约束即可求解t_y。如果I₁、I₂、I₄、I₅之中的任意两项为零，则

等于两个非零项平方和的倒数，否则

的值由下式9计算：

其中，

B＝(I₁I₅-I₂I₄)²；

假设t_y为正，由式8分别求出r₁₁、r₁₂、t_x、r₂₁和r₂₂的值，进而将r₁₁、r₁₂、t_x、r₂₁和r₂₂的值代入式6计算，可以得到x_c和y_c。

接下来可以根据旋转矩阵R的正交性继续求解R的其它元素，首先根据下式10计算r₁₃和r₂₃：

此处同样可以先假设r₁₃符号为正，而r₂₃的符号则与r₁₁r₂₁+r₁₂r₂₂的符号相反，旋转矩阵R第三行的元素分别为：

联立式3和式4可以得到下式12：

将上式12整理为矩阵形式，如下式13：

其中，E＝r₁₁X_w+r₁₂Y_w+t_x，F＝r₂₁X_w+r₂₂Y_w+t_y，G＝r₃₁X_w+r₃₂Y_w。同样采用最小二乘方法求解上式13的超定方程组即可获得有效焦距f、径向畸变系数k_i和平移分量t_z。

在本公开的一个可选实施方式中，需要说明的是，如上所述，在计算x_c和y_c时，首先假设t_y为正，由式8分别求出r₁₁、r₁₂、t_x、r₂₁和r₂₂的值，并将r₁₁、r₁₂、t_x、r₂₁和r₂₂的值代入式6计算得到x_c和y_c。如果x_c和y_c的符号与x_d和y_d的符号一致，则证明该假设是正确的。如果x_c和y_c的符号与x_d和y_d不一致，则需要将t_y的符号改为负，重新计算r₁₁、r₁₂、t_x、r₂₁和r₂₂的值，并将新的r₁₁、r₁₂、t_x、r₂₁和r₂₂的值代入式6重新计算x_c和y_c。

在本公开的一个可选实施方式中，需要说明的是，如上所述，采用最小二乘方法求解上式13的超定方程组时，如果求解得到的所述有效焦距f为负数，则需要回到式10将r₁₃和r₂₃改为相反的符号，并根据式11重新计算所述旋转矩阵R的第三行元素，进而更新式13的系数矩阵，即可获得正确的所述有效焦距f、所述径向畸变系数k_i和所述平移分量t_z。

至此，可以完成高速摄像机标定所需的全部内外参数的求解。可见基于共面点的分步线性标定方法的求解过程全部通过线性计算完成，标定过程简单快捷，并且在模型中考虑了三阶径向畸变，精度可以满足实验要求。

此外，在应用上述方法对高速摄像机进行标定的过程中，由于算法自身的限制，高速摄像机与二维标定板之间的关系还应满足以下条件：

(1)由于在中间变量的求解过程中平移分量t_y处在分母上(见式8)，因此在摆放二维标定板时应确保世界坐标系原点的实际成像位置远离成像平面的y轴，以保证t_y的取值不接近零。

(2)高速摄像机标定必须获得不同场深位置的图像信息，因此应保证标定点沿Z_c轴分布在不同位置，所以二维标定板不能垂直于摄像机光轴放置，即要求世界坐标系的Z_w轴不能平行于摄像机光轴，必须与光轴有一夹角。

(3)如果世界坐标系的原点在图像边缘或视野范围之外，那么在应用径向约束条件时将能够最大程度的消除透镜变形对有效焦距的影响，但这一点对于高速摄像机的标定并不是必须的。

在本公开的一个可选实施方式中，通过将两台高速摄像机呈光轴正交放置来实现在两台高速摄像机之间相互提供场深信息，即某一点在一台摄像机的画面中沿水平方向的位置即为该点在另一台摄像机中的场深。值得一提的是，应用该方法获取场深的前提是需要确保两台高速摄像机的光轴正交。

由式1至式5可知，观测点P的世界坐标(X_w,Y_w,Z_w)与其像点P_d的图像坐标(u,v)之间存在如下关系：

其中，矩阵

包含了从世界坐标系到图像坐标系的全部变换关系。

中所涉及的参数在上述高速摄像机标定过程中已经全部求解出来。由式18可知，若要实现从二维图像到三维空间的三维重构，即在已知(u,v)的情况下求解(X_w,Y_w,Z_w)，还需要获得z_c。z_c为观测点至摄像机光心的距离，即观测点的场深。

如图5所示，在正对高速摄像机的观察箱壁面上绘制有沿水平方向间隔为1cm的标尺，利用标尺可以对图像中某一点的水平位置进行估计，该估计值将作为三维重构中另一台高速摄像机的场深初值。在本实施方式中所用的世界坐标系的原点固定于观察箱的一角，各坐标轴的方向如图6所示。将1号摄像机拍摄的图像称为图像1，将2号摄像机拍摄的图像称为图像2。在图像1和图像2上获取的图像坐标分别用下标1和下标2表示，其余各坐标的表示方式同理。则对于图像1，根据图像坐标(u₁,v₁)以及在图像2上估计的场深初值

由式18即可求解得到图像1上的第一组世界坐标

对于图像2，根据图像坐标(u₂,v₂)以及在图像1上估计的场深初值

即可求解得到图像2上的第一组世界坐标

Y_w2和X_w1分别对应摄像机1和摄像机2的场深，因此令

并再次根据式24进行计算，即可求得图像1和图像2的第二组世界坐标

和

重复这一过程直到满足下式19：

在实际计算中一般迭代十余次即可获得满足式19的最终解。

在本公开的一个可选实施方式中，为了使实验结果更加精确，可以在对摄像机进行标定之前，先对两台摄像机的光轴的垂直度进行校准。在本实施方式中，可以将两台摄像机光轴正交转化为两台摄像机的光轴分别垂直于观察箱的对应壁面，并利用圆环标定板对其进行校准。

如图7所示，A、B、C和D点分别表示观察箱壁面上的四个点，O点表示摄像机光心，O₁点表示光轴与壁面的交点，A′、B′、C′、D′和O′点分别表示A、B、C、D和O点的成像，AD垂直于光轴，图中AB＝CD。首先考虑A点在倾斜前和倾斜后的成像，由针孔成像原理可获得下式20：

其中，A₁O₃表示A₁点到光轴的垂直距离。

则倾斜角θ对A点成像造成的影响可表示为下式21：

考虑AB倾斜后的成像，由针孔成像原理可得下式22：

其中，A₁′B₁′表示AB绕O₁点倾斜θ之后的成像，B₁′O₁′表示BO₁绕O₁点倾斜θ之后的成像，B₁O₂表示B₁点到光轴的垂直距离。

则AB倾斜后的成像可表示为下式23：

同理，CD倾斜后的成像可表示为下式24：

由式23和式24可得下式25：

Δ＝A′₁B′₁-C′₁D′₁∝θ 式25

当θ＝0时，Δ＝0，即如果壁面与光轴垂直，壁面上分别位于光轴两侧的两条相等线段的成像差为零。这是由于当光轴与物面存在夹角时，位于光轴两侧的物体放大倍数不同，对光轴垂直度的校准就利用了这条结论。在光轴垂直度的校准中选用如图8所示的圆环标定板作为标定物，将圆环标定板贴合观察箱壁面放置。上述结论只要求等长线段位于光轴两侧，故拍摄时只要保证光轴穿过圆环标定板的内圆即可，不一定要对准圆心。之后，在拍摄的圆环标定板的图像上测量AB、CD、EF和GH四条线段的长度。判定光轴与观察箱壁面的条件可表示为下式26：

其中，ε为误差阈值，因为印刷、测量等环节不可避免的存在误差，因此要求线段成像差完全为零是不可能的。在实际操作中只需低于某一阈值即可，阈值ε可以根据测量所需精度选取。如果Δ_N或Δ_H不能满足上式26，应根据式23和式24判断倾斜角θ并对摄像机进行相应调整，直到校准结果满足上述式26的判定条件。

接下来，可以开始进行图像坐标的提取和对应点匹配，以及三维数据重构。

如图9所示，本实施方式中，图像的分析处理通过Image Pro Plus6.0(IPP)完成。利用IPP的序列工具、测量工具、数据采集器以及宏录制功能，可以实现对蜻蜓飞行图像序列中特征点(观测点)的图像坐标的快速提取。在两台摄像机所拍摄的图像中分别提取相同特征点的图像坐标即完成了对应点的匹配，只有在两台摄像机中均有清晰成像的特征点才能够完成三维重构。

在本实施方式中，高速摄像机标定以及三维数据重构过程均已完成MATLAB(Matrix Laboratory)下程序的编写，可以通过程序计算结合IPP软件快速实现高速摄像机标定和三维数据重构。同时高速摄像机标定后如果不需要调整摄像机位置，此结果可以一直使用，无须重复标定。

下面结合附图详细说明本实施方式中所用两台高速摄像机的标定及结果。

如图10所示，对本实施方式中所用两台高速摄像机进行标定所借助的标定物为二维棋盘标定板。二维棋盘标定板的尺寸是600mmx900mm，棋盘图案中每一格尺寸为100mmx100mm。

根据上述高速摄像机与二维标定板之间应满足的所有条件，将二维棋盘标定板倾斜放置在观察箱中，并令尽量多的标定点落在视野范围内。

两台高速摄像机分别拍摄的标定图像如图11所示。提取图像中所有标定点的图像坐标(u,v)，应用上述方法对两台高速摄像机分别进行标定，获得摄像机内外参数初步标定结果，如下表1所示：

表1摄像机内外参数初步标定结果

需要注意的是，标定过程中世界坐标系建立在二维标定板上，表1中的旋转矩阵R和平移向量T相应为摄像机坐标系与标定板坐标系的变换矩阵。由于标定板是随意放置的，故该变换关系在之后的运算中没有任何实际意义。因为世界坐标系固定在观察箱的一角，所以需要再经过一次坐标变换，得到摄像机坐标系与观察箱坐标系的变换矩阵，两台高速摄像机的最终标定结果如下表2所示：

表2摄像机内外参数最终标定结果

基于两台高速摄像机的三维重构验证：

为了验证上述高速摄像机标定以及三维数据重构方法的准确性，利用二维棋盘标定板对其进行验证。将标定板倾斜放置在观察箱内，通过两台摄像机同时对其进行拍摄，得到图像如图12(a)和(b)所示。对比两幅图像，选取标定板中部在两幅图像中均有清晰成像的20个标定点为三维重构点。通过IPP对这20个标定点进行图像坐标提取，如图12(c)和(d)。三维重构的目标坐标系采用图8所示的观察箱坐标系。重构得到20个标定点在三维空间中的位置，如图13所示。棋盘标定板上相邻标定点间距为100mm，计算三维重构后相邻标定点的距离如表3所示，其中对标定点的编号方式同图12(c)和(d)。

表3三维重构的20个标定点的间距验证

P1	99.812	P2	100.275	P3	100.047	P4	100.562	P5
									101.134		100.538		100.993		100.514		100.930
P6	100.393	P7	99.975	P8	100.374	P9	100.028	P10
									100.487		100.780		100.204		100.387		100.454
P11	100.114	P12	99.903	P13	100.285	P14	100.426	P15
									101.107		100.828		100.913		100.876		100.671
P16	100.316	P17	100.465	P18	100.236	P19	100.180	P20

由表3可知，本实施方式对高速摄像机的标定以及三维数据重构方法具有较高的精度。重构后的标定点间距最大误差为1.134mm，即最大相对误差为1.13％。该误差包含了标定板印刷误差、摄像机标定误差、图像坐标提取误差和三维重构计算误差等，结合本实施方式中蜻蜓飞行观测实验的具体情况，认为该精度能够满足实验要求。

运用本公开技术方案进行实际拍摄的结果：

利用本公开的方法针对蜻蜓起飞过程进行了拍摄，获得了其翅膀和身体的运动参数，并基于此分析得出蜻蜓起飞过程的高升力机制，如图14和15所示。

在本公开中，所用的昆虫扑翼飞行实验观测系统可以方便的安置在实验室中，仅用两台高速摄像机对昆虫进行拍摄，实验成本相对较低，操作过程也较简单，对于活体昆虫没有特殊要求。本公开的基于两台高速摄像机的活体昆虫运动采集及数据重构方法采用仅考虑三阶径向畸变的摄像机模型，结合分布线性标定方法以及两台呈光轴正交的高速摄像机相互提供场深的方法，完成三维数据重构，操作过程简单明了，精度相对较高，并且可以完成昆虫长时间连续飞行测量。

本领域的技术人员应当理解，上述实施方式仅仅是为了清楚地说明本公开，而并非是对本公开的范围进行限定。对于所属领域的技术人员而言，在上述公开的基础上还可以做出其它变化或变型，并且这些变化或变型仍处于本公开的范围内。

Claims (8)

1.一种基于两台高速摄像机的活体昆虫运动采集及数据重构方法，其特征在于，两台所述高速摄像机的光轴垂直于供活体昆虫活动的观察箱的壁面，并且两台所述高速摄像机的光轴正交，

所述方法包括：

建立所述高速摄像机的坐标变换模型，其中，所述坐标变换模型包括：

刚体变换模型，用于从世界坐标系到摄像机坐标系的坐标变换；

透视投影模型，用于从摄像机坐标系到理想成像平面坐标系的坐标变换；

径向畸变模型，用于从理想成像平面坐标系到实际成像平面坐标系的坐标变换；

以及像素单位变换，用于从实际成像平面坐标系到图像坐标系的坐标变换；

通过建立所述坐标变换模型来确定需要求解的所述高速摄像机的参数；

通过基于共面点的分步线性标定方法求解所述高速摄像机的参数；

分别获取观测点在两台所述高速摄像机中的场深；以及

在两台所述高速摄像机所拍摄的图像中分别提取具有相同特征点的图像坐标来完成对应点的匹配，并根据所述高速摄像机参数和所述场深对所述观测点进行三维重构；

其中，通过将两台所述高速摄像机呈光轴正交放置，所述观测点在一台所述高速摄像机的图像中沿水平方向的距离是所述观测点在另一台摄像机中的场深；

所述方法还包括：对所述高速摄像机进行标定之前，将两台所述高速摄像机光轴正交转化为两台所述高速摄像机的光轴分别垂直于所述观察箱的对应壁面，并利用圆环标定板对两台所述高速摄像机的光轴的垂直度进行校准；

对所述高速摄像机进行标定所借助的标定物为二维棋盘标定板，摆放所述二维棋盘标定板时需确保世界坐标系原点的实际成像位置远离成像平面的y轴,且所述二维棋盘标定板相对于所述高速摄像机的光轴倾斜放置在所述观察箱中。

2.根据权利要求1所述的方法，其特征在于，

所述刚体变换模型由一个旋转矩阵R和一个平移向量T来描述，相应齐次坐标之间的关系如下式1和2：

其中，(x_c,y_c,z_c)表示所述观测点在所述高速摄像机的坐标系下的三维坐标，(X_w,Y_w,Z_w)表示所述观测点在世界坐标系下的三维坐标；

所述透视投影模型的相应齐次坐标之间的关系如下式3：

其中，f为所述高速摄像机的有效焦距，(x_u,y_u)表示所述观测点在成像平面坐标系下的理想成像坐标；

所述径向畸变模型是在所述透视投影模型的基础上加入畸变的影响，所述径向畸变模型仅考虑三阶径向畸变，相应齐次坐标之间的关系如下式4：

其中，(x_d,y_d)表示所述观测点在成像平面坐标系下的实际成像坐标，k_i表示各阶径向畸变系数，r表示像点到所述高速摄像机的光轴的距离；

所述像素单位变换是由物理长度向像素单位的变换，实际成像坐标与图像坐标之间的关系如下式5：

其中，(u,v)表示所述实际成像坐标在图像坐标系下的坐标，(u₀,v₀)表示所述成像平面坐标系的原点在图像坐标系下的坐标，dx表示单位图像元素在水平方向的尺寸大小，dy表示单位图像元素在竖直方向的尺寸大小，

表示图像元素比例因子。

3.根据权利要求2所述的方法，其特征在于，所述高速摄像机参数包括：

摄像机内部参数，包括所述成像平面坐标系的原点在图像坐标系下的坐标(u₀,v₀)、所述图像元素比例因子

所述有效焦距f和所述径向畸变系数k_i；以及

摄像机外部参数，包括所述旋转矩阵R和所述平移向量T。

4.根据权利要求3所述的方法，其特征在于，所述高速摄像机参数的求解方法包括：

将标定所用的世界坐标系OX_wY_wZ_w固定在二维标定板上，世界坐标系的X_w轴和Y_w轴平行于标定板，坐标系符合右手法则；

将所述二维标定板放置在所述高速摄像机的拍摄空间中，提取所述二维标定板的图像中各标定点的图像坐标(u,v)；

根据式5求解所述观测点在成像平面坐标系下的实际成像坐标(x_d,y_d)；

根据径向约束条件：

获得下式6：

将式6右端项的分子和分母同时乘以1/t_y，即将透镜变形对所述有效焦距f的影响考虑进来，获得下式7：

采用最小二乘方法求解式7的超定方程组并获得5个中间变量，所述中间变量分别为：

根据I₁、I₂、I₄、I₅求解t_y，如果I₁、I₂、I₄、I₅之中的任意两项为零，则

等于两个非零项平方和的倒数，否则

的值由下式9计算：

其中，

B＝(I₁I₅-I₂I₄)²；

假设t_y为正，由式8分别求出r₁₁、r₁₂、t_x、r₂₁和r₂₂的值，并将r₁₁、r₁₂、t_x、r₂₁和r₂₂的值代入式6计算得到x_c和y_c；

根据所述旋转矩阵R的正交性继续求解R的其它元素，首先根据下式10计算r₁₃和r₂₃：

此处同样先假设r₁₃符号为正，而r₂₃的符号则与r₁₁r₂₁+r₁₂r₂₂的符号相反，所述旋转矩阵R第三行的元素分别为：

联立式3和式4可以得到下式12：

将上式12整理为矩阵形式，如下式13：

其中，E＝r₁₁X_w+r₁₂Y_w+t_x，F＝r₂₁X_w+r₂₂Y_w+t_y，G＝r₃₁X_w+r₃₂Y_w，采用最小二乘方法求解上式13的超定方程组即可获得所述有效焦距f、所述径向畸变系数k_i和所述平移分量t_z。

5.根据权利要求4所述的方法，其特征在于，

假设t_y为正，由式8分别求出r₁₁、r₁₂、t_x、r₂₁和r₂₂的值，并将r₁₁、r₁₂、t_x、r₂₁和r₂₂的值代入式6计算得到x_c和y_c；如果x_c和y_c的符号与x_d和y_d的符号一致，则证明所述假设是正确的；

如果x_c和y_c的符号与x_d和y_d不一致，则将t_y的符号改为负，重新计算r₁₁、r₁₂、t_x、r₂₁和r₂₂的值，并将新的r₁₁、r₁₂、t_x、r₂₁和r₂₂的值代入式6重新计算x_c和y_c。

6.根据权利要求4所述的方法，其特征在于，如果求解得到的所述有效焦距f为负数，则回到式10将r₁₃和r₂₃改为相反的符号，并根据式11重新计算所述旋转矩阵R的第三行元素，进而更新式13的系数矩阵，即可获得正确的所述有效焦距f、所述径向畸变系数k_i和所述平移分量t_z。

7.根据权利要求4所述的方法，其特征在于，

所述二维棋盘标定板的尺寸为600mmx900mm，棋盘图案中每一格的尺寸为100mmx100mm。

8.根据权利要求1至7中任一项所述的方法，其特征在于，所述高速摄像机的拍摄区域设置在所述观察箱的中间区域。

Images