基于量子免疫克隆算法的三维人体运动跟踪方法

具体实施方式

本发明从图像序列的连续性出发，提出了一种基于量子免疫克隆算法的三维人体运动跟踪的新方法，即利用量子免疫克隆算法超强的搜索全局最优解的能力，进行三维的人体运动跟踪。本发明在对量子进化算法的一些改进的基础上，通过实验确立了相关的最优参数和统计特征，然后利用量子更新算子和量子交叉算子对状态参数进行诱导，并将跟踪结果和对应二维检测点进行了比较。详细描述如图1、图3和图4。

参照图1，本发明的具体实现过程如下：

步骤1，检测二维图像中人体关键关节点，并使用kalman滤波器对漏检点进行估计。

本发明使用的图像序列是从有人体运动存在的视频序列中获得的，总共660帧，由于本发明没有研究分割方法，本视频序列的人体运动没有自遮挡和被遮挡，而且由于对于人体躯干的旋转检测存在很大问题，本发明中的人体运动没有体现出人体躯干的旋转，同时为了减少人体模型维数，通常这里假设人体躯干上的关节点是没有自由度的。

从人的侧影信息中找到人体的关键关节点，是一项非常具有挑战性得工作，主要靠经验选取，无太多的理论，并且关键关节点的三维位置并不像运动捕获设备那样可以精确的获得。在检测关节点之前，首先做出以下假设：1)假设前两帧图像没有人体侧影自遮挡或者被遮挡，这个假设可以保证以后的有遮挡的关节点，可以根据一定规则预测出来；2)假设侧影的中心与三维骨架的根节点投影相差不远，这个假设，可以保证在任意情况下获取根节点。前景的人体侧影从图像序列分割出来，接着用数字形态学的方法，把人的侧影区域细化成单像素的图像骨架。其测过程如下：

(1.1)头部点检测：人体侧影中，头部的形状一般不会被遮挡，特征也最稳定，一般为圆形，因此，借鉴前人检测三维像素数据中头部的三维球壳模板，给出一个检测视频图像序列头部区域的2维同心圆模板。假设同心圆的内圆圆心为R₁，外圆圆心为R₂，同心圆模板的圆心沿着侧影的骨架S＝{s_i}搜索，i＝1，2，…，N，N为当前帧侧影区域骨架点的个数，计算落在内圆和外圆之间的轮廓C＝{c_j}的点数，j＝1，2，…，M为当前帧轮廓点的个数，当落入内外圆的点数最多时，s_i即为头部点，一般外圆半径要取内圆半径的2倍以上；

(1.2)根节点检测：假设人体侧影点集记作A＝{a_k}，点的个数为L，则中心点位置为根节点的经验检测法是取骨架点S＝{s_i}上的离人体侧影中心点最近的一点；

(1.3)肩关节与胯关节检测：假设头部点与根节点已经检测完毕，将左右肩关节与胯关节与整个躯干绑定，使其在检测过程中保持相对位置不变，把三维骨架的头部投影点与根节点投影点与图像上的对应点匹配好，让三维肩关节点投影到人的侧影区域；

(1.4)四肢端点检测：在侧影骨架线上，运用端点搜寻法，搜索骨架线的端点，一般说来，没有毛刺的情况下，通常很容易检测到四肢的位置，但是在噪声很大的情况下该方法将失效，于是给出一个禁搜区域，即在人头部点与根节点，两肩关节点之间划定一个矩形区域，该区域内搜索到的端点将被丢弃，这样就会极大地提高四肢点的检测准确率；

(1.5)膝关节与肘关节检测：首先将侧影骨架线在根节点上方的点集定义为上半身，其余的点定义为下半身，过下半部中心点作两脚的平行线，该平行线与骨架线的交点即为膝关节点，肘关节检测是寻找骨架线上的一点，使得该点到手和肩部的距离之差最小；

当人体运动发生自遮挡或者由于检测算法的不完备性导致某些关键点不能正确检测时，就必须实施遮挡点的预测机制，用经典卡尔曼滤波器进行遮挡点预测。

下面对卡尔曼滤波的基本理论进行简要描述，考虑离散时间线性随机动态系统：

x_k+1＝F_kx_k+Γ_kw_k

z_k＝H_kx_k+v_k

其中k∈N是时间指标，x_k∈Rⁿ为k时刻的系统状态向量，F_k是系统状态转移矩阵，w_k则为过程演化噪声，Γ_k为噪声矩阵，z_k∈R^m是k时刻对系统状态的量测向量，H_k为量测矩阵，v_k是量测噪声。

通过简单的外推就可以估计得到两帧之间缺失的关键关节点，同时对提取骨架和检测过程中的噪声进行了滤波，使得所有检测到的二维关键关节点的位置更合理。

步骤2，根据检测并预测处理后的二维人体关键关节点，建立一个虚拟的人体三维骨架模型。

(2.1)本发明中使用的人体骨架模型是用连接刚体连接的运动链，并使用一些先验的约束，以便生成身体的躯干部分，根结点和胯部，颈关节和肩部没有运动自由度，只是随着躯干运动；

(2.2)为了尽量降低搜索空间大小，通常采用的模型中，人体骨架共有21个自由度，包括3个整体的平移、6个胯关节、2个膝关节、6个肩关节、2个肘关节和2个颈关节；

(2.3)建立三维人体模型，首先根据人的解剖学比例大致确定人体各部分的比例，人体模型坐标系的中心为腹部节点，即为根节点，每个关节点都会有一个局部坐标系，其子节点坐标用矩阵M_P，C移动到父节点坐标：

M_P，C＝T(t_x，t_y，t_z)·R_z(γ)·R_y(β)·R_x(α)

其中T(t_x，t_y，t_z)表示子节点坐标沿x轴平移t_x、沿y轴平移t_y、沿z轴平移t_z，R_x(α)表示子节点坐标绕x轴旋转α度，R_y(β)表示子节点坐标绕y轴旋转β度，R_z(γ)表示子节点坐标绕z轴旋转γ度，它们的矩阵表示式为：

$T(t_x,t_y,t_z)=\left(\begin{array}{cccc}0& 0& 0& t_x\\ 0& 0& 0& t_y\\ 0& 0& 0& t_z\\ 0& 0& 0& 1\end{array}\right)$

$R_x\left(\alpha \right)=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \left(\alpha \right)& -\sin \left(\alpha \right)& 0\\ 0& \sin \left(\alpha \right)& \cos \left(\alpha \right)& 0\\ 0& 0& 0& 1\end{array}\right)$

$R_y\left(\beta \right)=\left(\begin{array}{cccc}\cos \left(\beta \right)& 0& \sin \left(\beta \right)& 0\\ 0& 1& 0& 0\\ -\sin \left(\beta \right)& 0& \cos \left(\beta \right)& 0\\ 0& 0& 0& 1\end{array}\right)$

$R_z\left(\gamma \right)=\left(\begin{array}{cccc}\cos \left(\gamma \right)& -\sin \left(\gamma \right)& 0& 0\\ \sin \left(\gamma \right)& \cos \left(\gamma \right)& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)$

对于子节点局部坐标系内的任意一点p＝(x，y，z，1)^T，在父节点局部坐标系内的位置为P＝(X，Y，Z，1)^T，它们之间的相互转化关系为：

$\left(\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right)=\left(\begin{array}{cccc}0& 0& 0& t_x\\ 0& 0& 0& t_y\\ 0& 0& 0& t_z\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{cccc}\cos \left(\gamma \right)& -\sin \left(\gamma \right)& 0& 0\\ \sin \left(\gamma \right)& \cos \left(\gamma \right)& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{cccc}\cos \left(\beta \right)& 0& \sin \left(\beta \right)& 0\\ 0& 1& 0& 0\\ -\sin \left(\beta \right)& 0& \cos \left(\beta \right)& 0\\ 0& 0& 0& 1\end{array}\right)$

$\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \left(\alpha \right)& -\sin \left(\alpha \right)& 0\\ 0& \sin \left(\alpha \right)& \cos \left(\alpha \right)& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\\ 1\end{array}\right)$

这样可以通过连续几个转化操作，就形成了一条运动链，表达出每个点的3D坐标与各个旋转角的关系，才能求出各点在世界坐标系中的位置，例如，左下肢各点坐标为：

P_LH＝T(t_x，t_y，t_z)R_0zR_0yR_0xT_LHR·[0，0，0，1]^T

P_LK＝T(t_x，t_y，t_z)R_0zR_0yR_0xT_LHRR_1zR_1yR_1xT_LKH[0，0，0，1]^T

P_LA＝T(t_x，t_y，t_z)R_0zR_0yR_0xT_LHRR_1zR_1yR_1xT_LKHR_2zR_2yR_2xT_LAK·[0，0，0，1]^T

P_LF＝T(t_x，t_y，t_z)R_0zR_0yR_0xT_LHRR_1zR_1yR_1xT_LKH

R_2zR_2yR_2xT_LAKR_3zR_3yR_3xT_LFK·[0，0，0，1]^T

其中，P_LH、P_LK、P_LA和P_LF分别代表左臀点，左膝点，左踝点以及左脚点，T(t_x，t_y，t_z)为根节点到世界坐标系原点的位移，R_ix为第i个节点绕x的旋转矩阵，R_iy为第i个节点绕y的旋转矩阵，R_iz为第i个节点绕z的旋转矩阵，下标0x、0y和0z代表躯干的整体旋转；T_LHR、T_LKH、T_LAK和T_LFA分别代表左臀部到根节点的初始平移、左膝盖到左臀节点的初始平移、左踝部到左膝盖的初始平移和左脚部到左踝部的初始平移；上肢点的求法类似于下肢点，这里不再一一列出，得到各关节点在世界坐标系内的坐标后，很容易根据各种投影模型推出各点在图像坐标系内的投影值，整体设计的三维人体模型，如图2所示。

步骤3，初始化种群，设置人体运动初始参数，并进行克隆操作。

(3.1)根据每个关键关节点的自由度，设置总体待估人体参数集X为：

X＝{x，y，z，θ₁，θ₂，…，θ₁₈}

其中x，y，z表示整体平移分别沿x、y、z轴的平移距离，θ₁，θ₂，…，θ₁₈表示待估计的人体关键关节点的旋转角度；由于人体运动的连续性，每两帧图像之间人体肢体变化角度一般不大，根据本图像序列相连帧的统计，设置两帧图像之间待估计的人体关键关节点的旋转角度θ₁，θ₂，…，θ₁₈的范围为[-15，+15]度；由于运动平移也很小，这里设置整体平移x，y，z的范围为[-15，+15]像素。

在估计过程中每一帧图像只需要得到两帧之间参数的差值ΔX：

ΔX＝{Δx，Δy，Δz，Δθ₁，Δθ₂，…，Δθ₁₈}

其中Δx，Δy，Δz表示两帧之间整体平移分别沿x、y、z轴平移的差值，而Δθ₁，Δθ₂，…，Δθ₁₈表示两帧之间待估计的三维人体状态角度的差值。

根据上述待估计的人体参数的设定，设置初始种群A(0)为：

其中每一列元素值是从+15到-15均匀赋值给这k+1个数字；每一行表示21个待估计的参数值；k表示每一个待估元素可以选取的状态值的大小；每一个初始元素对应设置有一个相等的概率1/k+1；

(3.2)进行克隆操作，生成克隆种群A′(k)。

根据初始化人体状态参数种群，为了扩张搜索空间，这里对初始种群进行克隆操作，表示为T_c^C：

$T_c^C\left(A\right)={\left(\begin{array}{cccc}{T}_c^C\left(a_1\right)& T_c^C\left(a_2\right),& ...& ,T_c^C\left(a_N\right)\end{array}\right)}^T$

其中，(i＝1，2，...，N)，N表示A(k)每一行的元素数，I_i为元素值为1的q_i维行向量，q_i为抗体a_i克隆后的规模，a_i表示A(k)中第i行的所有元素，q_i取为：

$q_i\left(k\right)=\mathrm{Int}(N_c\times \frac{D\left(a_i\left(k\right)\right)}{\underset{j=1}{\overset{N}{\Sigma }}D\left(a_j\left(k\right)\right)}),i=\mathrm{1,2},...,N$

Int(x)表示大于等于x的最小整数，而N_c是与克隆后的规模有关的设定值且满足N_c＞N，D(a_i(k))表示A(k)中每一个元素概率的大小，

克隆过后，初始种群根据每一个元素概率形成了多个搜索种群，这些种群总体定义为A′(k)：

A′＝{A，A′₁，A′₂，...，A′_N}

其中，A′_i表示经克隆操作后生成的新种群，N表示克隆后种群的规模。

步骤4，对进行克隆操作后的种群，进行免疫基因操作和免疫选择，生成新的种群。

(4.1)对克隆后的子群体进行免疫基因操作，生成免疫种群A″(k)；

经过克隆操作之后，人体状态参数的搜索空间根据人体关节点的距离亲和度扩张了数倍，为了防止人体状态函数的收敛速度慢和过早的陷于局部最优解，进行免疫基因操作，它主要包括克隆重组和克隆变异。使用量子更新算子实现克隆重组，使用量子交叉算子实现克隆变异；该量子更新算子，采用量子旋转门加速收敛，混沌变异操作，以防止早熟。

量子更新算子中的量子旋转门，表示为U(θ)：

$U\left(\theta \right)=\left(\begin{array}{cc}\cos \left(\theta \right)& -\sin \left(\theta \right)\\ \sin \left(\theta \right)& \cos \left(\theta \right)\end{array}\right)$

其中，θ为旋转更新的角度。

根据A′(k)中元素的概率，判断旋转的方向和角度，通过不断的更新和旋转最终实现更新操作。

量子更新算子中的混沌变异为了防止抗体朝着一个方面演化，而陷入局部最优，这里定义混沌变异，增加局部搜索能力。对于经上述量子旋转更新的抗体，以变异概率p_m随机选择一个或若干个位，实施如下操作：

${a^{\prime }}_{\mathrm{ji}}=\left(\begin{array}{ccc}{a}_{\mathrm{ji}}+(1-a_{\mathrm{ji}})\times {\mathrm{Logistic}}_{\mathrm{ji}}\left(k\right)& \mathrm{if}& \mathrm{rand}>0.5\\ a_{\mathrm{ji}}-a_{\mathrm{ji}}\times {\mathrm{Logistic}}_{\mathrm{ji}}\left(k\right)& \mathrm{if}& \mathrm{rand}\le 0.5\end{array}\right)$

其中a_ji为变异前第j个量子抗体中的第i个量子位的概率值，a′_ji为变异后第j个量子抗体中的第i个量子位的概率值。

混沌变异的变异尺度由两部分组成，即：

I.1-a_ji或a_ji能自动保证变异后抗体变量仍在[0，1]之间。

II.Logistic_ji(k)为Logistic映射的第k个序列值，其中Logistic映射的序列长度为q_i-1，即保留原抗体信息，对子群体内实施该变异操作，Logistic映射描述如下：

x_n+1＝μx_n(1-x_n)  (n＝0，1，2，...)

其中选取μ＝4，x_n＝a_ji。

(4.2)对种群A″(k)进行免疫选择，形成新的种群。

经过免疫基因操作之后，进行免疫选择，表示为T_S^C，它是从经克隆免疫基因操作后的各自子代和相应父代中选择优秀的抗体，从而形成新的种群。即：

$A(k+1)=T_s^C(A\left(k\right)\cup A^{\prime \prime }\left(k\right))$

根据亲合度，若有

b_i(k)＝{a″_ij(k)|maxD(a″_ij)j＝1，2，...，q_i-1}

使得

D(a_i(k))＜D(b_i(k))，i＝1，2，...，N

则

a_i(k+1)＝b_i(k)，i＝1，2，...，N

其中a″_ij表示经过更新操作后第i个量子抗体第j位的概率，D(a″_ij)表示a″_ij的亲和度，D(a_i(k))和D(b_i(k))表示a_i(k)的亲和度和b_i(k)亲和度，a_i(k+1)表示下一代的第i行元素。

即用b_i(k)替换A″(k)中的a_i(k)，如果不存在这样的b_i(k)，则a_i(k)不变，免疫选择完成后得到下一代抗体群A(k+1)。

步骤5，将种群中的状态参数代入三维人体模型，得到人体关键关节点三维坐标，利用二维关键关节点和三维投影点的距离，设计出距离函数。

在基于单目图像序列的人体运动估计中，在未知摄像机参数的情况下假设一个投影模型，在本发明中，使用了弱透视投影模型，预设一个根节点景深参数D，使三维人体模型能够投影到图像平面，与平面上检测的对应点相匹配。待估计的运动参数包括人体模型的整体平移，以及各关节旋转角度共21个参数X＝{x，y，z，θ₁，θ₂，…，θ₁₈}，人体三维模型上15个关键点P₁，P₂，…，P₁₅的三维坐标由这21个参数所决定，其中P_i＝(P_ix，P_iy，P_iz)投影到图像平面后的坐标为p_i＝(p_ix，p_iy)

$p_{\mathrm{ix}}=f\frac{P_{\mathrm{ix}}}{D}$

$p_{\mathrm{iy}}=f\frac{P_{\mathrm{iy}}}{D}$

检测到的图像内的关键点记为q_i＝(q_ix，q_iy)，令G(X)为二位检测点和三维投影点的对应加权和，则，

$G\left(X\right)=\underset{i=1}{\overset{17}{\Sigma }}{\left|\right|p_i-q_i\left|\right|}_2.$

步骤6，设计出相似度函数，计算距离亲和度，保存最佳解，恢复出三维人体姿态。

(6.1)建立了距离函数G(X)后，参数估计问题就是求一个非线性的优化问题：

X＝minG(X)

利用相似度函数X＝arg min G(X)，计算出二维关键关节点和三维投影点的距离的加权和，再求其最小值，每一代都存储最优解，本代中的最优解不大于全局最优解的不保存；设定一个极小值δ，如果全局最优解小于等于这个值，则终止计算，否则返回到步骤(3.2)。

(6.2)最后根据以前设计出的三维人体骨架模型，使用保存的最佳解代入模型，运用计算机图形学相关知识，对人体初始骨架模型进行旋转和平移，得到三维人体重建图。

本发明的效果可通过以下仿真进一步说明：

1.仿真内容：

采用本发明所提出的量子免疫克隆算法，进行三维人体姿态恢复实验。其中图像序列来自于自拍视频序列，采集到图像帧数660帧，图像序列间距1/24s。在估计过程中，通常每5帧采样一帧进行估计，得到132帧的结果。

首先用两个序列来测试头部点检测法，初始半径取得合适的情况下，头部检测准确率达到100％，超出了通常用来做初始点的要求。其他各关节点的检测基本符合本实验的要求。对于漏检点，通过Kalman滤波，通常可以很容易估计到检测不出来的特征点。

本发明将人体运动估计问题归结为确定性非线性优化问题，一般很难找到最优解，并且理论上，该问题会有2^N个局部最小解，N为旋转自由度的个数。初始值的选取以及运动约束是否合适将会极大影响优化的结果，在本实验中，人的初始姿态为身体站立，水平伸开双手，反映到数值表示为

X₀＝(167，96，0，

      15，0，0，-30，

      15，0，0，-30，

      0，0，

      0，0，0，10，

      0，0，0，-10)。

其中，前三个元素单位为像素，其他的均为角度，第17和21个元素为左右肘部的角度，分别设置为10度和-10度，即让胳膊微弯，第4和8个元素的设置保证大腿微前屈，第7与11个元素的设置保证小腿后伸。优化算法可以设置各参数分量的搜索范围，各参数的下边界记为LB，上边界为UB，总而言之，初值的设置是为了保证算法搜索过程中整个人的运动是合理的，不能出现肘部和膝部的反转。

根据亲和度评价状态值，存储最优解，恢复出人体三维姿态。

硬件平台为：Intel Core2 Duo CPU E6550@2.33GHZ、2GB RAM。软件平台为MATLAB 7.0。

2.仿真结果及分析

仿真结果如图3和图4所示。

图3给出部分估计结果，其中图3a是右腿前伸，上臂上举的动作；图3b是恢复出来的三维人体姿态从正面的视角观察；图3c是恢复出来的三维人体姿态从侧面的视角观察；图3d左腿前伸，上臂上举并回弯；图3e是恢复出来的三维人体姿态从正面的视角观察；图3f是恢复出来的三维人体姿态从侧面的视角观察；图3g是右腿前伸，双臂后伸的动作；图3h是恢复出来的三维人体姿态从正面的视角观察；图3i是恢复出来的三维人体姿态从侧面的视角观察。

从图3可以看出，利用量子免疫克隆算法可以估计得到满意的三维人体姿态。从图4可以看出，三维投影点和二维对应检测点的匹配很理想。并且所引入的算法运用到人体跟踪中计算代价较低，不容易陷入局部最优解。