GMM算法实现

西南交通大学，信息学院，[e-mail]资料整理1高斯混合模型（GMM）参数优化及实现龚勋（2010-11-13）1高斯混合模型概述高斯密度函数估计是一种参数化模型。有单高斯模型（SingleGaussianModel,SGM）和高斯混合模型（Gaussianmixturemodel，GMM）两类。类似于聚类，根据高斯概率密度函数（PDF，见公式1）参数的不同，每一个高斯模型可以看作一种类别，输入一个样本x，即可通过PDF计算其值，然后通过一个阈值来判断该样本是否属于高斯模型。很明显，SGM适合于仅有两类别问题的划分，而GMM由于具有多个模型，划分更为精细，适用于多类别的划分，可以应用于复杂对象建模。下面以视频前景分割应用场景为例，说明SGM与GMM在应用上的优劣比较：lSGM需要进行初始化，如在进行视频背景分割时，这意味着如果人体在前几帧就出现在摄像头前，人体将会被初始化为背景，而使模型无法使用；lSGM只能进行微小性渐变，而不可突变。如户外亮度随时间的渐变是可以适应的，如果在明亮的室内突然关灯，单高斯模型就会将整个室内全部判断为前景。又如，若在监控范围内开了一辆车，并在摄像头下开始停留。由于与模型无法匹配，车会一直被视为前景。当车过很长时间离去时，由于车停留点的亮度发生了很大的变化，因此已经无法与先前的背景模型相匹配；lSGM无法适应背景有多个状态，如窗帘，风吹的树叶。单高斯模型无法表示这种情况，而使得前背景检测混乱，而GMM能够很好地描述不同状态；l相对于单高斯模型的自适应变化，混合高斯模型的自适应变化要健壮的多。它能解决单高斯模型很多不能解决的问题。如无法解决同一样本点的多种状态，无法进行模型状态转化等。1.1单高斯模型多维高斯（正态）分布概率密度函数PDF定义如下：111(;,)exp[()()]2(2)TNp-∑=---∑∑xμxμxμ(1)注意与一维高斯分布不同，其中x是维数为d的样本向量（列向量），μ是模型期望，Σ是模型方差。对于单高斯模型，由于可以明确训练样本是否属于该高斯模型（如训练人脸肤色模型时，将人脸图像肤色部分分割出来，形成训练集），故μ通常由训练样本均值代替，Σ由样本方差代替。为了将高斯分布用于模式分类，假设训练样本属于类别C，那么，式(1)可以改为如下形式：111(/)exp[()()]2(2)TNCp-=---∑∑xxμxμ(2)式(2)表明样本属于类别C的概率大小。从而将任意测试样本ix输入式(2)，均可以得到一个标量(;,)iN∑xμ，然后根据阈值t来确定该样本是否属于该类别。²阈值t的确定：可以为经验值，也可以通过实验确定。另外也有一些策略可以参考：如文献[1]中采用搜索策略：首先令t＝0.7，以0.05为步长一直减到0.1左右，选择使样本变化最小的那个阈值做为最终t值，也就是意味着所选t值所构造的分类模型最稳定。²几何意义理解：根据单高斯分布PDF的含义我们可以知道，符合SGM分布的二维点在平面上应该近似椭圆形；相应地，三维点在空间中则近似于椭球状西南交通大学，信息学院，[e-mail]资料整理21.2高斯混合模型高斯混合模型是单一高斯机率密度函数的延伸，由于GMM能够平滑地近似任意形状的密度分布，因此近年来常被用在语音、图像识别等方面，得到不错的效果。例：有一批观察数据1{,...,}nX=xx，数据个数为n，在d维空间中的分布不是椭球状（如图1(a)），那么就不适合以一个单一的高斯密度函数来描述这些数据点的机率密度函数。此时我们采用一个变通方案，假设每个点均由一个单高斯分布生成（如图1(b)，具体参数jμ,jΣ未知），而这一批数据共由M（明确）个单高斯模型生成，具体某个数据ix属于哪个单高斯模型未知，且每个单高斯模型在混合模型中占的比例ja未知，将所有来自不同分布的数据点混在一起，该分布称为高斯混合分布。(a)(b)图1高斯混合模型图示，(a)表示所有样本数据；(b)表示已经明确了样本的分类从数学上讲，我们认为这些数据的概率分布密度函数可以通过加权函数表示：1()(;,)MijjijjjpNa==∑xxμΣ(3)上式即称为GMM，11Mjja==∑，其中()111;,exp()()2(2)TjjjjjjmjNp-⎡⎤=---⎢⎥⎣⎦xμΣxμΣxμΣ(4)表示第j个SGM的PDF。令(,,)jjjjqa=μΣ，GMM共有M个SGM模型，现在，我们就需要通过样本集X来估计GMM的所有参数：1(,,)TMqqΘ=⋅⋅⋅。样本X的概率公式为：()()1111()log;,log;,NMjjjjjiNMjjjjijlXNNaa====Θ==∑∏∑∑xμΣxμΣ(5)西南交通大学，信息学院，[e-mail]资料整理32采用EM估计GMM的参数通常采用EM算法（期望值最大，ExpectationMaximum）根据(5)式对GMM参数进行估计。具体推导过程在很多文献可以参阅，这里不详述，我们主要介绍采用EM估计GMM参数的方法。描述估计高斯混合分布的EM算法（参考“!!EM算法.ppt”）：²初始值：方案1：协方差矩阵0jΣ设为单位矩阵；每个模型比例的先验概率：01/jMa=；均值0jμ设为随机数；方案2：由k均值（k-mean）聚类算法对样本进行聚类，利用各类的均值作为0jμ，并计算0jΣ，0ja取各类样本占样本总数的比例。²算法流程：（1）估计步骤（E-step）：令ja的后验概率为：(;)(;)(;)jjijjiMllilNENabaaΘ=Θ=Θ∑xxx，1in≤≤，1jM≤≤(6)注意：写代码实现公式(6)时，对于每个SGM分别用公式(4)计算每个样本点ix在该模型下的概率密度值(;)jiNΘx，对于所有样本，得到一个n*1的向量，计算M次，得到*nM的矩阵，每一列为所有点在该模型下的概率密度值（PDF）；实现(;)MllilNaΘ∑x时，需要针对每个点计算在各个SGM的概率值总和。公式(6)可以通过如下matlab代码描述：functionprob=GaussPDF(Data,Mu,Sigma)[dim,N]=size(Data);Data=Data'-repmat(Mu',N,1);prob=sum((Data*inv(Sigma)).*Data,2);prob=exp(-0.5*prob)/sqrt((2*pi)^dim*(abs(det(Sigma))+realmin));（2）最大化步骤（M-step）：Ø更新权值：'1NijijNba==∑(7)Ø更新均值：'11nijiijnijibb===∑∑xμ(8)Ø更新方差矩阵：'''11()()nTijijijijnijibmmb==--=∑∑xxΣ(9)注意，有些文献把方差矩阵当作对角阵处理，即认为2()jdiags=Σ，故采用如下公式估计2s：西南交通大学，信息学院，[e-mail]资料整理4''211()()nTijijijinijibmmsb==--=∑∑xx(10)经本人实验验证，公式(10)会给最终结果带来误差，故在运算速度能保证的情况下建议采用公式(9)估计方差矩阵。（3）收敛条件方案1：不断地迭代E和M步骤，重复更新上面的三个值，直到'()()lXlXeΘ-Θ，通常510e-=，()lXΘ通过公式(5)计算，'()lXΘ表示更新参数后计算的值；方案2：不断地迭代E和M步骤，重复更新上面的三个值，直到参数的变化不显著，即'eΘ-Θ，'Θ为更新后的参数，通常510e-=。注：实验表明方案1和方案2效果接近，但方案2明显运算量要小，故推荐使用。3源码3.1单高斯模型下面代码实现了SGM，并实现了人脸肤色检测。其中图像处理、矩阵运算采用了openCV库函数/*****************************************************************************SingleGaussianModelforskincolorextractionParam:img--inputimagetoextractthefaceregionskinImg--result*****************************************************************************/voidCSkinColor::RunSGM(IplImage*img,IplImage**skinImg){if(img==NULL)return-1;////////////////////////////////////////////////////////////////////////////以下参数一组(117.4361,156.5599)来自源码light2，与文章《王航宇:基于YCbCr高斯肤色模型的//人脸检测技术研究》相同，另一组来自源码“肤色检测正式版”（103.0056,140.1309）doubleM[]={103.0056,140.1309}/*{117.4361,156.5599}*/;//M为肤色在YCbCr颜色空间的样本均值(Cb,Cr)，经验值doubleC[2][2]={{160.1301,12.1430},//C为肤色相似度模型的协方差矩阵，同上为经验值{12.1430,299.4574}};//注：因为运算仅需要该矩阵的逆矩阵值，故该值没有使用，仅作参考doubleinvC[2][2]={0.0077,-0.0041,-0.0041,0.0047};//Ct为C的逆矩阵值，由matlab计算而得//////////////////////////////////////////////////////////////////////////IplImage*pImg=img;doubleCrMean=0,CbMean=0,YMean=0;//1颜色转换：BGR-YCrCbIplImage*imgYCrCb=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U,3);//YCrCb图像西南交通大学，信息学院，[e-mail]资料整理5cvCvtColor(pImg,imgYCrCb,CV_BGR2YCrCb);//第0,1,2层分别为Y,Cr,CbIplImage*imgY=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U,1);//YCrCb图像IplImage*imgCr=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U,1);//YCrCb图像IplImage*imgCb=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U,1);//YCrCb图像IplImage*imgY32=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_32F,1);//YCrCb图像IplImage*imgCr32=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_32F,1);//YCrCb图像IplImage*imgCb32=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_32F,1);//YCrCb图像cvSplit(imgYCrCb,imgY,imgCr,imgCb,NULL);cvConvert(imgY,imgY32);cvConvert(imgCr,imgCr32);cvConvert(imgCb,imgCb32);////////////////////////////////////////////////////////////////////////////2根据SigleGaussianModel计算颜色模型IplImage*PCbCr=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_