0 引言

人脸特征点定位是人脸分析过程的重要步骤，定位信息能为人脸图像深入处理与分析提供重要前提依据，在人脸识别^[1]、人脸跟踪^[2]、表情分类^[3]等领域已存在广泛应用且应用意义较高。

近些年来，回归技术在模式识别领域成功应用为本文研究工作提供基础理论条件：将人脸特征点定位看成回归问题，直接对人脸图像外观进行投影学习回归函数，从而产生目标输出^[4]。传统回归方法并未强调形状约束条件，只保证初始形状符合人脸形状约束条件。在复杂人脸图像中，眼睛、眉毛和嘴巴部分因为各种实际表情(挑眉、大笑、眯眼等)，或是遮挡(头发、眼镜等)会造成特征点定位误差并超出人脸合理区域范围，因此本文选择添加姿态估计来执行形状约束操作。针对整个人脸执行姿态估计只能判断人脸是否产生旋转，无法获取每个动作导致的五官具体变化的问题，本文针对人脸特征点进行分块操作，对各个部分姿态分别进行估计，再利用估计好的姿态进行后续定位。

目前，姿态估计方法包括基于模型的方法^[12]和基于学习的方法^[13]。利用几何模型来表示物体结构，通过提取特征在模型和图像之间建立对应关系从而实现姿态估计。这种方法需要匹配多种局部特征，而受光照强度、噪声大小等影响无法准确地提取局部特征。基于学习的方法从事先标定的不同姿态训练样本中学习对应关系，将学习所得回归函数应用于样本，所得结果作为对应样本姿态估计。该方法采用全局特征，具有较好鲁棒性，然而其只能将姿态划分到预先定义姿态类别中，不能保证估计精度和连续性。

为了解决上述问题，本文提出应用特征姿态核回归模型的姿态估计算法。首先，通过训练能够学习到输入变量(图像特征)与输出变量(姿态)的联合概率分布函数(匹配函数)；然后，通过求解联合概率分布函数(匹配函数)在输出空间的最大值，获取目标姿态估计；同时改进级联回归算法，将姿态估计作为形状约束，并对形状特征采样添加先验信息，进行特征点精确定位。

1 核回归模型

回归分析是确定两种或两种以上变量之间相互依赖定量关系的统计分析方法，包括逻辑回归、线性回归、非线性回归、非参数回归等。主要内容是从一组数据出发，确定某些变量之间的定量关系表达式，这种表达式称为回归方程，利用所求回归方程进行预测和控制，如式(1)所示：

${y_i} = m({\mathit{\boldsymbol{x}}_i}) + {\varepsilon _i}{\kern 1pt} $

(1)

其中:m(x)是回归函数，变量y是向量，变量x是向量，i范围为1~N(N代表特征维数)，ε是随机误差。核回归函数作为模型，对于给定训练数据集{(x_i, y_i)}_i=1^N，m(x)可表示为一系列核函数在训练数据集或其子集上的线性组合，如式(2)所示：

$m(\mathit{\boldsymbol{x}}) = E(\left. \mathit{\boldsymbol{y}} \right|\mathit{\boldsymbol{x}}){\rm{ = }}\sum\limits_{i = 1}^N {{\alpha _i}{k_i}\left( {\mathit{\boldsymbol{x}}{\kern 1pt} ,{\mathit{\boldsymbol{x}}_i}} \right)} $

(2)

核函数k_i(x, x_i)表示输入变量x之间关系，α_i通过线性方程求解。

2 特征姿态核回归模型及求解 2.1 特征姿态核回归模型

以COFW(Caltech Occluded Faces in the Wild)人脸数据库为实例，如图 1所示，总共包括29个特征点，共分为7块区域(左眉、右眉、左眼、右眼、鼻子、嘴巴、下巴)，最易受到干扰部分是眉毛、眼睛、嘴巴。其中，左眼以(9、11、13、14) 4个点构成四边形表示，17号点则用(13、14)两点连线中点替换。以最上方点为一号点逆时针排序。左眉、右眉、右眼部分类似，嘴巴部分则是用(23、24、25、28)表示，(26、27)两点在(25、28)两点连线上，与(25、28)两点距离固定。

本文模型输入变量是包含目标部分的人脸图像数据，输出变量是四边形中1、2、3号点对应角度以及1、2号点连线和2、3号点连线的长度，以此来确定四边形形状及姿态。传统模型的输出变量视为一维变量，对于一维输出变量，一般采用一个回归函数分别估计每一维输出。在学习1号点角度回归函数时，同一个输出值可能会对应多幅图像，而这些图像的2、3号点角度不同，训练过程无法反映图像姿态完整信息。

为解决算法时间复杂度问题，重新定义训练数据集{(x_i, y_i)}_i=1^N。其中，y_i是训练图像(或图像特征)x_i对应姿态向量。定义核函数k_Y(y, y_j)表示输出向量y之间的关系，将式(2) $\sum\limits_{i = 1}^N {{\alpha _i}{k_i}\left( {\mathit{\boldsymbol{x}}{\kern 1pt} ,{\mathit{\boldsymbol{x}}_i}} \right)} $执行加权操作得到式(3)：

$M(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{y}}) = \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {{\beta _{ij}}{k_X}\left( {\mathit{\boldsymbol{x}}{\kern 1pt} ,{\mathit{\boldsymbol{x}}_i}} \right){k_Y}\left( {\mathit{\boldsymbol{y}}{\bf{,}}{\mathit{\boldsymbol{y}}_j}} \right)} } $

(3)

其中：k_X(x, x_i)和k_Y(y, y_j)分别定义输入与输出空间核函数，表示输入变量之间和输出变量之间关系；β_ij通过训练获取；M(x, y)是匹配函数，表示对于某个输入x、输出为y的可能性。当定义M(x, y)取值在时，可将其视为x和y的联合概率分布函数，式(3)可改成式(4)所示：

$\begin{array}{l} M(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{y}}) = \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {{\beta _{ij}}{k_X}\left( {\mathit{\boldsymbol{x}}{\kern 1pt} ,{\mathit{\boldsymbol{x}}_i}} \right){k_Y}\left( {\mathit{\boldsymbol{y}}{\kern 1pt} ,{\mathit{\boldsymbol{y}}_j}} \right) = } } \\ \quad \sum\limits_{i = 1}^N {\left[ {{\beta _{i1}},{\beta _{i2}}, \cdots ,{\beta _{iN}}} \right]} \left[ {\begin{array}{*{20}{c}} {{k_Y}\left( {\mathit{\boldsymbol{y}},{\mathit{\boldsymbol{y}}_1}} \right)}\\ {{k_Y}\left( {\mathit{\boldsymbol{y}},{\mathit{\boldsymbol{y}}_2}} \right)}\\ \vdots \\ {{k_Y}\left( {\mathit{\boldsymbol{y}},{\mathit{\boldsymbol{y}}_N}} \right)} \end{array}} \right]{k_X}\left( {\mathit{\boldsymbol{x}}{\kern 1pt} ,{\mathit{\boldsymbol{x}}_i}} \right) = \\ \quad \sum\limits_{i = 1}^N {{\lambda _i}} {k_X}\left( {\mathit{\boldsymbol{x}}{\kern 1pt} ,{\mathit{\boldsymbol{x}}_i}} \right) \end{array}$

(4)

定义φ(y)=[k_Y(y, y₁), k_Y(y, y₂), …, k_Y(y, y_N)]^T，A=(β_ij)∈R^N×N，则λ_i=β_i·φ(y)表示输出空间的权重。本文模型通过在输出空间定义核函数k_Y(y, y_j)加权到传统核回归模型核函数上，可同时表示输入变量和输出变量关系，通过一次学习就能够获取整个模型回归函数，不需要针对每一维输出进行单独训练来提高训练效率。

2.2 图像特征及核函数选择

为有效提高速度和识别率，本文选取图像特征为像素差特征(Pixel-Difference Feature, PDF)^[11]。每个像素灰度差特征h_{p1, p2}(I)就是定义在图像中p₁、p₂位置上的两个像素点差值，赋值方法如式(5) :

${h_{{p_1},{p_2}}}(I) = I({p_1}) - I({p_2})$

(5)

其中：I(p)表示图像I在位置p的灰度值。像素差特征采用线性插值法，依据目标形状而不是原始图像坐标对像素点进行索引，如图 2所示。这样可确保特征几何不变性而简化回归模型，在训练过程中收敛速度更快，在处理充足训练样本时，计算过程依然简便且鲁棒性好。

输入和输出空间核函数均采用高斯核函数形式表达，具体如式(6)所示:

$\left\{ \begin{array}{l} {k_X}\left( {\mathit{\boldsymbol{x}}{\kern 1pt} ,{\mathit{\boldsymbol{x}}_i}} \right) = \exp \left( {\frac{{ - \left\| {\mathit{\boldsymbol{x}}{\kern 1pt} - {\mathit{\boldsymbol{x}}_i}} \right\|_X^2}}{{2{\sigma _X}}}} \right)\\ {k_{{\kern 1pt} Y}}\left( {\mathit{\boldsymbol{y}}{\kern 1pt} ,{\mathit{\boldsymbol{y}}_j}} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \,{\kern 1pt} {\kern 1pt} {\kern 1pt} \exp \left( {\frac{{ - \left\| {\mathit{\boldsymbol{y}}{\kern 1pt} - {\mathit{\boldsymbol{y}}_j}} \right\|_Y^2}}{{2{\sigma _Y}}}} \right) \end{array} \right.$

(6)

其中: ‖·‖_X和‖·‖_Y是分别定义在输入和输出空间的第二范数。x为通过训练确定的每幅图像的像素值差，y为姿态向量，即输出3个角度值和2个距离值。核函数参数设定如式(7)所示:

$\left\{ \begin{array}{l} {\sigma _X} = 0.1 \times \mathop {\max }\limits_i \left\| {\mathit{\boldsymbol{x}}{\kern 1pt} - {\mathit{\boldsymbol{x}}_i}} \right\|_X^2\\ {\sigma _{{\kern 1pt} Y}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = 0.1 \times \mathop {\max }\limits_i \left\| {\mathit{\boldsymbol{y}}{\kern 1pt} - {\mathit{\boldsymbol{y}}_j}} \right\|_{{\kern 1pt} Y}^2 \end{array} \right.$

(7)

2.3 模型求解

定义全部训练数据匹配函数矩阵M∈R^N×N元素为：

${M_{ij}} = M\left( {{\mathit{\boldsymbol{x}}_i},{\mathit{\boldsymbol{y}}_j}} \right) = \left\{ {\begin{array}{*{20}{c}} 1\\ 0 \end{array}} \right.\begin{array}{*{20}{c}} ,\\ , \end{array}\begin{array}{*{20}{c}} {i = j}\\ {i \ne j} \end{array}$

(8)

$\mathit{\boldsymbol{M}} = {\mathit{\boldsymbol{K}}_X} \cdot \mathit{\boldsymbol{A}} \cdot {\mathit{\boldsymbol{K}}_Y}$

(9)

其中：M为N阶单位矩阵，K_X、K_Y均为N阶方阵。通过训练得到系数矩阵如式(10)所示：

$\mathit{\boldsymbol{A}} = \mathit{\boldsymbol{K}}_X^{ - 1} \cdot \mathit{\boldsymbol{M}} \cdot \mathit{\boldsymbol{K}}_Y^{ - 1}$

(10)

定义ψ(x)=[k_Y(x, x₁), k_Y(x, x₂), …, k_Y(x, x_N)]^T，则匹配函数最终表达为式(11)所示：

$\mathit{\boldsymbol{M}}{\bf{(}}\mathit{\boldsymbol{x}},\mathit{\boldsymbol{y}}) = \mathit{\boldsymbol{\psi }}\left( \mathit{\boldsymbol{x}} \right) \cdot \mathit{\boldsymbol{A}} \cdot \mathit{\boldsymbol{\varphi }}(\mathit{\boldsymbol{y}})$

(11)

对于未知图像${\mathit{\boldsymbol{\hat x}}}$，姿态估计如式(12)：

$\mathit{\boldsymbol{\hat y}} = \arg \max m\left( {\mathit{\boldsymbol{\hat x}},\mathit{\boldsymbol{y}}} \right)$

(12)

3 利用姿态估计的级联回归 3.1 传统级联回归

令S=[x₁, y₁, …, x_L, y_L]^T是被拉伸为列向量、在人脸图像I中的L个特征点的(x, y)坐标，即人脸形状。令S^t是人脸形状S的当前估计(S⁰表示初始人脸形状)，S_i表示人脸图像数据库中对人脸形状基准标定(面部特征点实际位置)，S表示训练样本集中S_i的平均值。

在回归前，首先对训练样本的人脸形状进行预处理，将当前输入形状S^t投影到S，即对S^t执行相似变换，从而保证回归过程中人脸形状保持不变。相似变换通过最小化两者之间绝对距离来实现：

$\left| {\mathit{\boldsymbol{\bar S}} - (\mathit{\boldsymbol{s}}_i^t\mathit{\boldsymbol{R}}_i^t{\mathit{\boldsymbol{S}}^t} + \mathit{\boldsymbol{T}}_i^t)} \right|$

(13)

式(13)中：s_i^t、R_i^t、T_i^t分别表示所求相似变换尺度变换矩阵、旋转变换矩阵、平移变换矩阵。

获取标准化形状之后，采用递进式回归策略，训练T个有序树结构相关弱回归器R=(R¹, R², …, R^T)。

输入I，S

for t=1 to T

  x^t=h^t(S^t-1, I)    //计算形状索引特征

  δS=R^t(x^t)    //估计回归量

  S^t=S^t-1+δS    //更新估计S

End

输出最终估计S^T

在回归阶段中，树状回归器R^t的每个叶子节点输出是形状增量ΔS^t=R^t(h^t(I, S^t-1))，再把所得形状增量ΔS^t与上次迭代所得形状估计S^t-1相加得到更新的新形状S^t作为当前形状估计。给定粗略初始形状估计S⁰，针对人脸形状估计S^T依照上述方式通过迭代方法求解：

${\mathit{\boldsymbol{S}}^T} = {\mathit{\boldsymbol{S}}^0} + \sum\limits_{t = 1}^T {\Delta {\mathit{\boldsymbol{S}}^t}} $

(14)

训练每个回归器R^t时，目标是使损失函数S_i^t=S_i^t-1+R^t(x_i^t)达到最小值。

${R^t} = \mathop {\arg \min }\limits_R \sum\limits_i {d(R({x_i}),{\mathit{\boldsymbol{S}}_i})} $

(15)

在训练好R^t后，利用计算得到S_i^t用于回归器训练。如果回归器不能再减小误差，则停止训练^[10]。

3.2 利用姿态估计的改进级联回归

针对级联回归方法进行改进，利用姿态估计作为形状约束。对每块区域进行回归时，因为姿态已知(即：四边形形状已确定)，只需确定任意两个点位置便可完成每个部分整体定位。初始位置设定是所有训练数据平均值，因为不需要计算所有点位置，可以更快降低估计误差从而收敛并减少级联个数，提高计算效率。

考虑人脸具有结构性，像素点间距离越近，几何意义也更近。例如眼角与眉毛附近的像素灰度差的判别性要强于眉毛与下颌附近像素灰度差，如图 3所示。

为了确保得到形状索引特征距离相近，在采样过程中补充指数先验进行距离约束。

${P_t}(\mathit{\boldsymbol{u}},\mathit{\boldsymbol{v}}) \propto {{\rm{e}}^{ - \alpha \left\| {\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{v}}} \right\| + \beta t}}$

(16)

式(16)中：‖u-v‖是像素点对之间距离，t是当前回归器数量，α、β是用来控制距离的参数。初始回归时，只需进行粗略定位，采样像素点对距离较大；随着回归器增加，为了对特征点进行精确定位，选择距离将不断缩小，从而可以应对人脸表情细微变化。

4 实验结果与分析

为验证本文方法对于遮挡和形变目标特征点定位的准确性，选择权威人脸数据库进行测试，包括LFPW(Labeled Face Parts in the Wild)、HELEN、COFW(Caltech Occluded Faces in the Wild)。本文选择2000张图片作为训练样本，其余330张图片作为测试样本；COFW包含1007张姿态表情变化较大、遮挡较为严重的图像，与LFPW具有相同基准标定。本文采用845张LFPW图像和500张COFW图像作为训练样本，其他507张COFW人脸图像作为测试样本。实验条件是Windows 7、i5-4440 3.1 GHz、8 GB内存配置的PC机上使用Matlab 2014b对比分析本文算法性能。

4.1 核函数选择实验

本文评价指标包括两个：1)平均误差，即估计值与真值误差绝对值的平均值；2)正确率，即估计值与真值误差绝对值小于某一阈值的百分比。

为评价输入和输出变量的核函数的选取对本文方法的性能影响，选取高斯核函数和Sigmoid核函数进行比较，实验在一维数据上进行，即只对单一角度结果进行计算。对比实验结果如表 1所示，对于输入变量，采用高斯核函数的姿态估计结果优于Sigmoid核函数；对于输出变量，在平均误差方面，采用高斯核函数较好，在2°和5°范围内的误差百分比则是Sigmoid核函数占优。结合两种核函数本身性质，输入变量间差别较大，采用高斯核函数进行表示较为合适，而输出变量变化范围有限，采用Sigmoid核函数进行表示，可提高整体估计的准确性，但相较高斯核函数其估计精度有所下降。因此，本文选择输入和输出空间均采用高斯核函数。

4.2 图像特征选择实验

实验在二维数据上进行，即对两种角度结果进行计算，实验结果如表 2所示, 其中：傅里叶描述子(Fourier Descriptor, FD)是一种基于频域变换的形状表示算法; 方向梯度直方图(Histogram of Oriented Gradients, HOG)是一种在计算机视觉和图像处理中用来进行物体检测的特征描述子；局部二值模式(Local Binary Pattern, LBP)特征是一种用来描述图像局部纹理特征的算子。

实验结果表明采用不同图像特征对于姿态估计结果有较大影响。在进行人脸部位姿态估计时，选择对于目标姿态变化较为敏感的图像特征。

4.3 特征点定位实验

为验证本文方法对存在遮挡以及形变的人脸特征点定位性能，采取评价指标是平均误差(Error)、失败率(Failure)、测试速度(帧率)。其中，定义最终估计特征点位置与数据库中标定特征点间的距离与每幅图片中双瞳的距离之比为误差，而误差大于0.1，则会被判定为失败，从而计算失败率。

将特征点划分为7个分块(3.1节已述)，每个分块回归器在训练过程中采用固定参数设置：第一层回归器数量为10，第二层回归器数量为300。每个回归器选择形状索引特征数量及所构建弱回归器层数是5。选择显示形状回归算法(Explicit Shape Regression, ESR)^[5]、鲁棒姿势级联回归算法(Robust Cascade Postural Regression, RCPR)^[6]、监督下降法(Supervised Descent Method, SDM)^[7]、局部二值特征(Local Binary Feature, LBF)算法^[8]及级联高斯过程回归树(cascade Gaussian Process Regression Tree, cGPRT)算法^[9]作为本文比较算法。

从表 3可见，本文算法在LFPW和HELEN数据库上与对比方法接近。在COFW数据库中，针对存在遮挡及形变较大图像，相较于其他算法定位精度和成功率有较大提高。

由于本文算法主要目标是处理存在遮挡的图像，LFPW与HELEN中存在遮挡图像较少，因此在这两个数据库中的实验主要用来验证分块回归和初始粗定位性能。下述将介绍在COFW数据库上的对比实验，选择ESR^[5]、RCPR^[6]及本文算法进行分块而不添加形状约束作为比较算法。

针对图像进行分块操作，本文算法相较于ESR算法在定位精度方面有较大提高，但仍落后于RCPR算法。加入姿态估计作为形状约束后的本文算法定位精度有提高，在回归未开始时就已经明显领先于其他算法并最终输出最准确结果，失败率降至14%左右，能够处理存在遮挡及形变人脸图像。本文算法处理速度已达到280 frame/s水平，能够满足实时处理标准要求。

5 结语

本文应用姿态估计分块特征点定位算法，通过获取图像特征与目标姿态间的匹配函数，并求取匹配函数最大值来获取目标姿态估计效果。为了提高算法对于复杂人脸图像和遮挡的定位精度，对人脸特征点进行分块操作，利用姿态估计作为每块形状约束，对于部分位置点进行回归，降低回归器数量来提高计算效率，并引入形状索引特征采样先验，使算法能够高效采样到判别性强的特征来构建回归器，从而提高算法实时性能。实验表明选择合适图像特征，本文算法可提供精确姿态估计，从而获得优于现有算法定位成功率和准确度，具有较强鲁棒性以及较高计算速度，能够达到实时处理水平。