在沟通交往中，人的长相、表情以及头部姿态提供了丰富的信息. 对人脸眉毛、眼角和嘴唇等关键点的定位是很多应用领域关键的一部分，比如头部姿态预测、人脸跟踪、人脸识别^[1-3]、人脸表情识别^[4]和医疗诊断的面部分析^[5]等. 基于二维图像的技术对光照、姿态和表情变化敏感，且易于伪造，具有安全性隐患. 随着三维采集技术和设备的发展，基于三维人脸的技术得到了更多的关注和研究.

统计形状模型是计算机视觉领域一项常用的技术，可以用于定位二维图像中人脸、肢体和器官等目标物体上预定义的特征点. Cootes等提出基于二维图像的主动形状模型^[6]、主动外观模型^[7]和约束局部模型^[8]（constrained local model, CLM）. 这些方法在对局部表面建模、训练和匹配优化策略部分使用不同的方法. 一些学者尝试使用基于二维图像的方法解决三维人脸特征点定位问题，Morency等^[9]基于CLM模型提出CLM-Z框架，该方法结合二维图像的灰度信息和三维人脸的深度信息，使用正则化特征点均值漂移算法^[10]作为模型匹配算法. Cheng等^[11]提出结合三维人脸深度、天顶角和方位角信息的CLM模型，在人脸深度图上定位特征点. 这些方法仅仅使用了人脸的深度信息，无法很好地处理姿态变化.

基于三维人脸点云的统计形状模型方法通常需要提取一个关键点集. Nair等^[12]建立三维点分布模型（point distribution model, PDM），该方法结合形状指数和曲度指数提取出特征点，生成特征点组合，通过最小化与模型形状的残差来实现定位. Perakis等^[13]提出大的姿态变化下检测特征点的统计人脸特征点模型，使用形状指数和旋转图提取出关键点，通过普氏分析和旋转图相似度筛选出特征点. Creusot等^[14]提出基于三维人脸的机器学习方法，计算训练集人脸每个特征点的10个局部描述符，为每个特征点训练一个分类器. 融合分类器的响应检测出关键点，使用随机抽样一致算法（random sample consensus, RANSAC）进行特征点选择. Sukno等^[15]提出可以处理人脸缺失遮挡的方法，为每个特征点统计出按旋转图相似度从高到低排列，包含有效候选点需要提取的候选点数量. 提取出数量足够的候选点，使用基于RANSAC的方法定位特征点. 这些方法需要对关键点集中特征点的候选点集进行嵌套遍历，计算开销随着特征点数的增加迅速增长.

本文提出基于特征融合约束局部模型的三维人脸特征点定位算法. 为了更好地利用三维人脸网格的局部信息，该算法结合人脸的深度和表面形状信息，为每个特征点分别训练一个深度信息分类器和形状信息分类器. 使用训练得到的深度图人脸检测器的特征点平均位置作为算法的初始化位置，利用正则化特征点均值漂移算法进行模型拟合，实现特征点定位. 该方法由于使用模型拟合代替对候选点的嵌套遍历，与基于关键点检测的定位方法相比，在待定位较多特征点时的时间开销更小.

约束局部模型分为以下3部分：形状模型、局部模型和模型拟合. 形状模型和局部模型分别用来描述物体的结构形状和局部外观信息；模型拟合利用局部模型获得的响应，在形状模型的约束下调整参数向量，使当前模型拟合到最优位置.

大多数的形变模型匹配方法对非刚性物体的形变采用线性优化策略，本文使用点分布模型对非刚性形变进行线性建模，并与全局的刚性形变结合. 将特征点的坐标连接成一个向量，表示二维图像中的形状：

(1) ${{\varOmega }} = [{{{\omega }}_1},{{{\omega }}_2},\cdots,{{{\omega }}_n}]. $

式中： ${{{\omega }}_i} = [{x_i},{y_i}]$ 表示第 $i$ 个特征点的位置坐标， $n$ 为待定位特征点数. 使用普氏分析^[16]对训练集中形状进行对齐，对对齐后的形状集使用主成分分析（principal component analysis, PCA），可以得到形状模型：

(2) ${{\varOmega }}({{p}}) = s{{R}}({\bar{ \varOmega }} + {{\varPhi q}}) + {{t}}. $

式中： ${\bar{ \varOmega }}$ 表示库集形状的普氏平均； ${{\varPhi }}$ 为特征向量组成的正交基矩阵； ${{p}} = [s,{{R}},{{q}},{{t}}]$ 是PDM的参数，包括比例因子 $s$、旋转矩阵 ${{R}}$、非刚性形变向量 ${{q}}$ 和平移向量 ${{t}}$. 使用高斯分布来近似参数 ${{p}}$ 的分布，形状先验概率如下：

(3) $P({{p}}) \propto N({{q}};0,{{\varLambda }}). $

式中： ${{\varLambda }} = {\rm{diag}}\; [{\lambda _1},\cdots,{\lambda _m}] $ 为PCA的特征值 ${\lambda _i}$ 组成的对角矩阵.

局部模型由 $n$ 个局部分类器组成，在待检测特征点周围区域获取局部外观信息，对局部外观信息进行建模得到响应函数，计算特征点在图像像素点对齐的概率作为局部模型在该点的响应：

(4) ${{{R}} _i}({{{\omega }}_i};{{{I}}}) = P({l_i} = 1|{{{\omega }}_i},{ I}). $

式中： ${l_i} \in {\rm{\{ 1}}, - {\rm{1\} }}$ 是一个离散随机变量，表示第 $i$ 个特征点是否对齐到正确的位置； ${R_i}$ 表示第 $i$ 个局部分类器的归一化响应函数； $P({l_i} = 1|{{\bf{\omega }}_i},{ {{I}}})$ 表示图像 ${ {{I}}}$ 中第 $i$ 特征点在位置 ${ \omega _i}$ 匹配的概率.

CLM拟合的目标是搜索PDM参数 $p$，使得所有特征点的对齐误差最小化，正则化形式如下：

(5) $\mathop {\arg \min }\limits_{{p}} \;[L({{p}}) + \sum\limits_{i = 1}^n {{E_i}({{{\omega }}_i};{ I})} ]. $

式中： $L({{p}})$ 为形变的惩罚项， ${E_i}({{{\omega }}_i};{ I})$ 为第 $i$ 个特征点在图像 ${ I}$ 位置 ${{{\omega }}_i}$ 的对齐误差. 式（5）可以解释为在所有特征点正确对齐的条件下，最大化 p 的概率. 假设特征点分类器之间相互独立，则有下式：

(6) $P({{p}}|\{ {l_i} = 1\} _{i = 1}^n,{ I}) \propto P({{p}})\prod\limits_{i = 1}^n {P({l_i} = 1|{{{\omega }}_i},{ I})} . $

式（5）中的正则项和对齐误差函数可以表示为

(7) $L({{p}}) = - \ln\; \{ P({{p}})\}, $

(8) ${E_i}({{{\omega }}_{{i}}};{ I}) = - \ln \;\{ P({l_i} = 1|{{{\omega }}_{{i}}},{ I})\} . $

CLM是局部方法，依赖于初始参数估计. 假设PDM参数的初始估计为 ${{{p}}_0}$，需要找到参数更新值 $\Delta {{p}}$，使得 ${{{p}}_0} + \Delta {{p}}$ 接近最优解 ${{{p}}_{\rm{*}}}$. 使用正则化特征点均值漂移法作为模型拟合方法，式（5）的最小二乘解可以表示为

(9) $\mathop {\arg \min }\limits_{\Delta {{p}}} \;(||{{{p}}_0} + \Delta {{p}}||_{{{{\varLambda}} ^{ - 1}}}^2 + ||{{J}}\Delta {{p}} - {{v}}|{|^2}). $

式中： ${{J}}$ 为特征点位置 ${{\varOmega }}$ 关于参数 ${{p}}$ 的雅克比矩阵， ${{v}} = [{{{v}}_1},\cdots,{{{v}}_n}]'$ 为基于局部响应的均值漂移向量. 用 $\{ {{\varPsi} _i}\} $ 表示特征点候选区域，对候选位置的计算约束在一个局部矩形区域， ${{{\varPsi}} _i}$ 表示区域内所有像素位置，对局部响应使用高斯核密度估计，则均值漂移向量表示为

(10) ${{{v}}_i} = \sum\limits_{{\eta _i} \in { \varPsi _i}} {\frac{{{\pi _{{\eta _i}}}N({{\omega }}_i^{\rm{c}};{\eta _i},\rho {{I}})}}{{\sum\nolimits_{{\xi _i} \in { \varPsi _i}} {{\pi _{{\xi _i}}}N({{\omega }}_i^{\rm{c}};{\xi _i},\rho {{I}})} }} - {{\omega }}_i^{\rm{c}}} . $

式中： ${\pi _{{\eta _i}}} = p({l_i} = 1|{\eta _i},{ I})$， ${{\omega }}_i^{\rm{c}}$ 表示第 $i$ 个特征点的当前预测位置，参数 $\rho $ 可以按文献[17]的方法从训练集中推测出来.

使用吉洪诺夫高斯-牛顿正则化方法，导出 $\Delta {{p}}$ 的解，使用 $\Delta {{p}}$ 对当前预测位置 ${{{\varOmega }}^{\rm{c}}}$ 进行更新，参数的更新如下：

(11) $\Delta {{p}} = - {({{{J}}^{{{\rm{T}}}}}{{J}} + \rho {{{\varLambda }}^{ - 1}})^{ - 1}}(\rho {{{\varLambda }}^{ - {{1}}}}{{p}} - {{{J}}^{{{\rm{T}}}}}{{v}}), $

(12) ${{\varOmega }} \approx {{{\varOmega }}^{\rm{c}}} + {{J}}\Delta {{p}}. $

如图1所示， 提出融合深度信息和网格局部信息的局部模型. 对三维人脸点云进行均匀重采样，得到分辨率为1 mm的三维点云. 将点云转化为深度图以及网格化，在三维网格上的点与深度图的人脸像素点之间建立一一对应关系. 分别使用深度信息和三维网格信息训练分类器，融合2个分类器响应作为局部模型的响应.

深度信息表示三维人脸到镜头的距离信息，参考文献[9]的方法，使用logistic回归表示深度信息响应：

(13) ${r_\lambda } = [{{1 + \exp \;({{\alpha }}{C_\lambda }({{\omega }},{{ I}_{{\rm{depth}}}}) + \beta )}}]^{-1}. $

式中： ${{\alpha }}$ 为logistic回归器的回归系数； $\beta $ 为截距； ${C_\lambda }$ 为特征点 $\lambda $ 的局部分类器的输出，使用支持向量机作为局部分类器， ${C_\lambda }$ 可以表示为

(14) ${C_\lambda }({{\omega }},{{ I}_{{\rm{depth}}}}) = {{{w}}_\lambda }V({{\omega }},{{ I}_{{\rm{depth}}}}) + {b_\lambda },$

其中 ${{{w}}_\lambda }$ 和 ${b_\lambda }$ 为权重和阈值， $V({{\omega }},{{ I}_{{\rm{depth}}}})$ 为以 ${{\omega}} $ 为中心的 $11 \times 11$ 的图像块的灰度连成的归一化向量.

参考文献[14]，计算人脸网格的10个局部描述符. 包括6个曲率相关的标量描述符：最大曲率 ${k_1}$、最小曲率 ${k_2}$、高斯曲率 $K$、平均曲率 $H$、形状指数 $S{\rm{_I}}$ 和对数均方差曲率 $L{\rm{_C}}$，如下：

(15) $K = {k_1}{k_2},$

(16) $H = \frac{{{k_1} + {k_2}}}{2},$

(17) $S _{\text{I}} = (2/\pi) \;\arctan \left[ \left({k_1} + {k_2}\right)/\left({k_1} - {k_2}\right) \right],$

(18) $L_{\text{C}} = (2/\pi) \;\log \sqrt {\left( {k_1^2 + k_2^2} \right)/2} .$

2个与局部平面相关的标量描述符：局部平面距离和局部体积，表示网格点到其周围点形成的平面的最近距离以及该点与平面构成立方体的体积. 2个直方图描述符：旋转直方图^[18]和球面直方图，分别统计一个点周围柱面和球面空间的点分布情况.

对于训练中的每个特征点，标量描述符值的分布可以用参数化的类条件概率密度分布近似. 使用2种概率密度函数（probability density function, PDF），对形状指数描述使用逆高斯分布函数，其他描述符使用高斯分布. 对于网格上的点 ${{\omega }}$ 以及描述符-特征点对 $(d,\lambda )$，可以通过用该点描述符的原始值 ${x_d}({{\omega }})$ 计算得分. 描述符-特征点得分使用下式计算：

(19) $s_\lambda ^d({{\omega }}) = \frac{{{\rm{pdf}}_\lambda ^d({x_d}({{\omega }}))}}{{\mathop {\max }\limits_{{x_d}} \;({\rm{pdf}}_\lambda ^d({x_d}))}}. $

式中： ${\rm{pdf}}_\lambda ^d$ 为描述符 $d$ 在特征点 $\lambda $ 处的概率密度； ${x_d}({{\omega }})$ 为描述符 $d$ 在点 ${{\omega }}$ 处的值； $s_\lambda ^d({{\omega }})$ 为 $d$ 在 ${x_d}({{\omega }})$ 处的概率密度与概率密度函数最大值的比值，表示特征点 $\lambda $ 在该点对齐的相对概率.

对于直方图描述符得分，需要计算描述符与平均直方图的差值并映射到一个标量上. 对于点 ${{\omega}} $，定义该点邻域为人脸网格上距离该点<5 mm的区域，非邻域为距离为15~45 mm的区域，将映射问题转化为LDA二分类问题：

(20) ${x_d}({{\omega }}) = {{w}}_\lambda ^{\rm{T}}({{{x}}'_d}({{\omega }}) - {{\bar{ x}}'_{d,\lambda }}). $

式中： ${{\bar{ x}}'_{d,\lambda }}$ 为特征点 $\lambda $ 直方图描述符的平均特征向量； ${{{w}}_\lambda }$ 为LDA训练出的分类器权重，根据式（19）对 ${x_d}({{\omega }})$ 进行归一化得到直方图描述符得分.

将10个描述符的得分级联成向量 ${{{s}}_\lambda }$，利用LDA训练邻域和非邻域2个点集的最佳分割，并将邻域和非邻域的响应值归一化到[0, 1]内，得到权重系数 ${{{u}}_\lambda }$，融合得到每个特征点形状信息分类器的响应：

(21) ${S_\lambda } = {{{s}}_\lambda }{{{u}}_\lambda }. $

人脸深度信息是三维点云在某个方向上的投影，仅仅使用深度信息无法有效地将人脸其他特征点和周围点区分开来. 三维人脸的表面具有更丰富的形状信息，不随着姿态变化发生形变，提取的特征对表情变化更加鲁棒. 本文融合深度信息和形状信息，使用 $R_\lambda ^1$ 表示特征点 $\lambda $ 在深度图某个点对齐的概率，为该点在深度信息响应 ${r_\lambda }$ 与 ${S_\lambda }$ 的平均值：

(22) $R_\lambda ^1 = ({r_\lambda } + {S_\lambda })/2. $

对1.3节的模型拟合算法进行修改. 使用深度图人脸检测器进行人脸检测，将训练集中特征点相对于人脸检测器的平均位置作为特征点初始化位置. 由于深度图内、外眼角距离大约为30像素，使用眼角距离一半的15像素作为矩形搜索区域的大小. 模型拟合算法的流程如下.

1）使用深度图人脸检测器初始化特征点位置，根据式（2）对模型参数 ${{P}}$ 进行初始化.

2）根据式（13）、（21）为每个特征点分别计算搜索窗口内的 ${r_\lambda }$ 和 ${S_\lambda }$，按式（22）对响应进行融合得到 ${R_\lambda }$.

3）根据式（10）计算均值漂移向量，按式（11）计算PDM参数更新 $\Delta {{P}}$，更新参数 ${{P}}$.

4）当 $\Delta {{P}} > {\rm{MIN\_STEP}}$ 且 $i < {\rm{MAX\_ITER}}$ 时，转到步骤2)，否则算法结束. 其中MIN_STEP为收敛条件，MAX_ITER为最大迭代次数，本文中MAX_ITER取10.

如图2所示，定位人脸的眼角、鼻尖、鼻梁、鼻翼、鼻下、嘴角、唇上和唇下，共14个点. 为了评估特征点检测方法的效果，使用2个人脸数据库：FRGC v2.0^[19]和Bosphorus^[20]. FRGC v2.0数据库是大型的公共人脸数据库，包含466个对象，共4 007张点云，人脸样本有着接近正面的姿态和不同的表情. Bosphorus数据库是另一个较新的三维人脸数据库，包含105个对象（60个男人、45个女人），共4 666张人脸点云，包含各种姿态表情和不同程度的遮挡. 每个人脸样本包含一个二维彩色图像、三维点云和标记好的24个人脸特征点.

根据文献[13]的方法，从FRGC v2.0中选取300张正面人脸作为训练集，手动标记出人脸的14个特征点，为每个特征点分别训练深度信息分类器和形状信息分类器. 从FRGC v2.0中的149个对象里选出975张正面人脸，作为测试集.

使用平均误差、误差标准差和检测成功率，对特征点定位的效果进行量化评估. 特征点的平均误差表示测试集算法定位的特征点与真实特征点位置的距离平均值，误差标准差为训练集中这个特征点误差距离的标准差. 个体检测误差是一个样本14个特征点误差的平均值. 总体检测误差是训练集所有个体检测误差分布的均值和标准差. 总体检测成功率是训练集中个体检测误差在某个阈值中的比率，本文选取阈值为6 mm.

如表1所示，本文方法的总体平均误差和标准差分别为3.12和1.37 mm，总体检测成功率为97.3%. 其中，每个特征点的平均误差为2.48~4.12 mm. 本文方法对内眼角点的定位效果最佳，检测成功率分别是97.0%和96.5%；对下巴点的定位效果最差，检测成功率只有76.3%.

为了验证特征融合的有效性，对仅基于深度信息和形状信息的CLM分别开展实验. 仅基于深度信息的CLM算法的总体误差和标准差为5.07和3.69 mm，仅基于局部形状信息的CLM算法的总体误差和标准差为3.76和2.34 mm. 可以看出，融合深度信息和形状信息后的特征点效果大大提高.

为了评估算法在姿态和表情变化下的效果，开展以下2个实验.

1）将FRGC v2.0的测试集按照表情强度分为中性、温和和极端3个子集进行实验，计算3种表情强度下子集的总体平均误差和检测成功率. 3个子集的总体检测成功率分别为97.6%、97.4%和95.5%. 可以看出，在表情强度变化下，总体检测成功率没有发生很大的变化，该方法在表情变化下具有较好的鲁棒性.

2）使用Bosphorus库中不同姿态的人脸，对该方法在姿态变化下的鲁棒性进行测试. Bosphorus库包含7个尺度的偏航旋转、4个尺度的俯仰旋转、2个角度的交叉旋转. 随机抽取100张姿态为轻微上仰的人脸作为数据集DB1，100张姿态为轻微俯视的人脸作为数据集DB2，100张偏航角为10°的人脸作为数据集DB3. DB1、DB2和DB3的检测成功率分别为94%、95%和89%，相较于FRGC v2.0测试集97.3%的检测成功率没有很大的下滑，可得该方法在小幅度的姿态变化下具有较好的鲁棒性.

1）与其他算法的精度对比.

Creusot等^[14]将每个三维网格分类器的响应融合，提取出人脸关键点，使用RANSAC算法搜索与模型匹配的候选点组合. 文献[14]算法的总体平均误差为4.95，误差标准差为3.33. Sukno等^[15]将人脸点集按照旋转图相似度从高到低进行排列，统计出训练集中对于每个特征点要提取有效候选点，需要从高到低遍历的点的数量为50. 提取出每个特征点的候选点，使用RANSAC算法实现特征点定位. 文献[15]算法的总体平均误差为3.44. Camgöz 等^[21]采用基于深度图的特征点定位方法. 该算法使用监督下降方法，训练一个通用下降方向序列和偏差项序列，更新特征点位置，实现特征点定位. 该方法仅仅用三维人脸的深度信息定位人脸的12个特征点，总体平均误差为4.97 mm，其中鼻翼点的平均检测误差分别为7.77和7.61 mm. 实验结果表明，在FRGC v2.0人脸集下，本文的特征点平均误差和总体平均误差都小于其他方法.

2）与其他算法的效率对比.

Creusot等^[14-15]在定位特征点都使用RANSAC算法，需要反复选择数据中的随机子集，找到与模型匹配最佳的组合. Creusot等^[14]使用阈值筛选出一组稀疏关键点集，算法的效率与关键点集的稀疏程度成正比，与定点精度成反比. 本文对FRGC v2.0数据库的点云进行均匀重采样，预处理后的每个点云包含大约4万个点，算法在i5处理器的Windows操作系统上使用C++进行编码，平均运行时间为1.05 s，使用GPU加速算法，效率可以提高几十倍. 文献[14]方法将FRGC v2.0数据库的点云稀释到平均包含3 000多个点，平均运行时间为1.18 s. 文献[15]方法给出的定位14个特征点的算法时间为31.5 s，该方法将点云稀释到平均包含4万个点.

（1）提出基于三维网格、融合三维人脸深度信息和局部形状信息的约束局部模型. 分别为每个特征点训练深度信息分类器和形状信息分类器，有效利用三维人脸丰富的局部表面信息，将特征点与其他点区分开来.

（2）正则化特征点均值漂移算法使用对初始形状进行迭代回归的方法，代替对候选特征点穷举遍历，生成候选特征点组合，避免了嵌套循环造成的随着待定位特征点数量迅速增长的时间开销.

本文基于FRGC v2.0数据库和Bosphorus数据库，对该方法进行实验评估. 实验结果表明，本文方法在精度上取得了较好的效果，对于表情变化和小幅度的姿态变化具有较好的鲁棒性.