高斯混合模型(GMM)参数优化及实现

西南交通大学，信息学院，xgong@swjtu.cn 资料整理 1 高斯混合模型（GMM）参数优化及实现 龚 勋 （2010-11-13） 1 高斯混合模型概述 高斯密度 函 关于工期滞后的函关于工程严重滞后的函关于工程进度滞后的回复函关于征求同志党风廉政意见的函关于征求廉洁自律情况的复函 数估计是一种参数化模型。有单高斯模型（Single Gaussian Model, SGM）和高斯混合模型 （Gaussian mixture model，GMM）两类。类似于聚类，根据高斯概率密度函数（PDF，见公式 1）参数的 不同，每一个高斯模型可以看作一种类别，输入一个样本 x，即可通过 PDF计算其值，然后通过一个阈值 来判断该样本是否属于高斯模型。很明显，SGM 适合于仅有两类别问题的划分，而 GMM 由于具有多个 模型，划分更为精细，适用于多类别的划分，可以应用于复杂对象建模。 下面以视频前景分割应用场景为例，说明 SGM与 GMM在应用上的优劣比较： l SGM 需要进行初始化，如在进行视频背景分割时，这意味着如果人体在前几帧就出现在摄像头前， 人体将会被初始化为背景，而使模型无法使用； l SGM 只能进行微小性渐变，而不可突变。如户外亮度随时间的渐变是可以适应的，如果在明亮的室 内突然关灯，单高斯模型就会将整个室内全部判断为前景。又如，若在监控范围内开了一辆车，并在 摄像头下开始停留。由于与模型无法匹配，车会一直被视为前景。当车过很长时间离去时，由于车停 留点的亮度发生了很大的变化，因此已经无法与先前的背景模型相匹配； l SGM 无法适应背景有多个状态，如窗帘，风吹的树叶。单高斯模型无法 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示这种情况，而使得前背 景检测混乱，而 GMM能够很好地描述不同状态； l 相对于单高斯模型的自适应变化，混合高斯模型的自适应变化要健壮的多。它能解决单高斯模型很多 不能解决的问题。如无法解决同一样本点的多种状态，无法进行模型状态转化等。 1.1 单高斯模型 多维高斯（正态）分布概率密度函数 PDF定义如下： 11 1( ; , ) exp[ ( ) ( )] 2(2 ) TN p - å = - - - å åx μ x μ x μ (1) 注意与一维高斯分布不同，其中 x是维数为d 的样本向量（列向量）， μ是模型期望， Σ是模型方差。 对于单高斯模型，由于可以明确训练样本是否属于该高斯模型（如训练人脸肤色模型时，将人脸图像 肤色部分分割出来，形成训练集），故 μ通常由训练样本均值代替， Σ由样本方差代替。为了将高斯分布 用于模式分类，假设训练样本属于类别C，那么，式(1)可以改为如下形式： 11 1( / ) exp[ ( ) ( )] 2(2 ) TN C p - = - - - å åx x μ x μ (2) 式(2)表明样本属于类别C的概率大小。从而将任意测试样本 ix 输入式(2)，均可以得到一个标量 ( ; , )iN åx μ ， 然后根据阈值 t来确定该样本是否属于该类别。 ² 阈值 t的确定：可以为经验值，也可以通过实验确定。另外也有一些策略可以参考：如文献[1]中采用 搜索策略：首先令 t＝0.7，以 0.05为步长一直减到 0.1左右，选择使样本变化最小的那个阈值做为最 终 t值，也就是意味着所选 t值所构造的分类模型最稳定。 ² 几何意义理解：根据单高斯分布 PDF的含义我们可以知道，符合 SGM分布的二维点在平面上应该近 似椭圆形；相应地，三维点在空间中则近似于椭球状 西南交通大学，信息学院，xgong@swjtu.cn 资料整理 2 1.2 高斯混合模型 高斯混合模型是单一高斯机率密度函数的延伸，由于 GMM 能够平滑地近似任意形状的密度分布， 因此近年来常被用在语音、图像识别等方面，得到不错的效果。 例：有一批观察数据 1{ ,..., }nX = x x ，数据个数为 n，在 d 维空间中的分布不是椭球状（如图 1(a)），那 么就不适合以一个单一的高斯密度函数来描述这些数据点的机率密度函数。此时我们采用一个变通 方案 气瓶 现场处置方案 .pdf气瓶 现场处置方案 .doc见习基地管理方案.doc关于群访事件的化解方案建筑工地扬尘治理专项方案下载 ， 假设每个点均由一个单高斯分布生成（如图 1(b)，具体参数 jμ , jΣ 未知），而这一批数据共由M （明确） 个单高斯模型生成，具体某个数据 ix 属于哪个单高斯模型未知，且每个单高斯模型在混合模型中占的比例 ja 未知，将所有来自不同分布的数据点混在一起，该分布称为高斯混合分布。 (a) (b) 图 1 高斯混合模型图示，(a)表示所有样本数据；(b)表示已经明确了样本的分类 从数学上讲，我们认为这些数据的概率分布密度函数可以通过加权函数表示： 1 ( ) ( ; , ) M i j j i j j j p Na = = åx x μ Σ (3) 上式即称为 GMM， 1 1 M j j a = =å ，其中 ( ) 11 1; , exp ( ) ( )2(2 ) T j j j j j jm j N p -é ù= - - -ê úë û x μ Σ x μ Σ x μ Σ (4) 表示第 j个 SGM的 PDF。 令 ( , , )j j j jq a= μ Σ ，GMM共有M 个 SGM模型，现在，我们就需要通过样本集 X 来估计 GMM的所有 参数： 1( , , ) T Mq qQ = ××× 。样本 X的概率公式为： ( ) ( ) 11 1 1 ( ) log ; , log ; , N M j j j j ji N M j j j j i j l X N N a a == = = Q = = åÕ å å x μ Σ x μ Σ (5) 西南交通大学，信息学院，xgong@swjtu.cn 资料整理 3 2 采用 EM估计 GMM的参数 通常采用 EM算法（期望值最大，Expectation Maximum）根据(5)式对 GMM参数进行估计。具体推导 过程在很多文献可以参阅，这里不详述，我们主要介绍采用 EM估计 GMM参数的方法。 描述估计高斯混合分布的 EM算法（参考“!!EM算法. ppt 关于艾滋病ppt课件精益管理ppt下载地图下载ppt可编辑假如ppt教学课件下载triz基础知识ppt ”）： ² 初始值： 方案 1：协方差矩阵 0jΣ 设为单位矩阵；每个模型比例的先验概率： 0 1/j Ma = ；均值 0jμ 设为随机数； 方案 2：由 k 均值（k-mean）聚类算法对样本进行聚类，利用各类的均值作为 0jμ ，并计算 0jΣ ， 0ja 取各 类样本占样本总数的比例。 ² 算法流程： （1）估计步骤（E-step）： 令 ja 的后验概率为： ( ; ) ( ; ) ( ; ) j j i j j i M l l i l N E N a b a a Q = Q = Qå x x x ， 1 i n£ £ ，1 j M£ £ (6) 注意：写代码实现公式(6)时，对于每个 SGM分别用公式(4)计算每个样本点 ix 在该模型下的概率密度 值 ( ; )j iN Qx ，对于所有样本，得到一个 n *1 的向量，计算M 次，得到 *n M 的矩阵，每一列为所有点在该 模型下的概率密度值（PDF）；实现 ( ; ) M l l i l Na Qå x 时，需要针对每个点计算在各个 SGM 的概率值总和。公 式(6)可以通过如下 matlab代码描述： function prob = 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）： Ø 更新权值： ' 1 N iji j N b a == å (7) Ø 更新均值： ' 1 1 n ij ii j n iji b b = = = å å x μ (8) Ø 更新方差矩阵： ' ' ' 1 1 ( )( )n Tij i j i ji j n iji b m m b = = - - = å å x x Σ (9) 注意，有些文献把方差矩阵当作对角阵处理，即认为 2( )j diag s=Σ ，故采用如下公式估计 2s ： 西南交通大学，信息学院，xgong@swjtu.cn 资料整理 4 ' ' 2 1 1 ( ) ( )n Tij i j i ji n iji b m m s b = = - - = å å x x (10) 经本人实验验证，公式(10)会给最终结果带来误差，故在运算速度能保证的情况下建议采用公式(9)估 计方差矩阵。 （3）收敛条件 方案 1：不断地迭代 E和M步骤，重复更新上面的三个值，直到 '( ) ( )l X l X eQ - Q < ，通常 510e -= ， ( )l X Q 通过公式(5)计算， '( )l X Q 表示更新参数后计算的值； 方案 2：不断地迭代 E和 M步骤，重复更新上面的三个值，直到参数的变化不显著，即 ' eQ - Q < ， 'Q 为 更新后的参数，通常 510e -= 。 注：实验表明方案 1和方案 2效果接近，但方案 2明显运算量要小，故推荐使用。 3 源码 3.1 单高斯模型 下面代码实现了 SGM，并实现了人脸肤色检测。其中图像处理、矩阵运算采用了 openCV库函数 /***************************************************************************** Single Gaussian Model for skin color extraction Param: img -- input image to extract the face region skinImg -- result *****************************************************************************/ void CSkinColor::RunSGM(IplImage *img, IplImage **skinImg) { if (img == NULL) return -1; ////////////////////////////////////////////////////////////////////////// // 以下参数一组(117.4361,156.5599)来自源码 light2，与文章《王航宇:基于 YCbCr 高斯肤色模型的 // 人脸检测技术研究》相同，另一组来自源码“肤色检测正式版”（103.0056, 140.1309） double M[]={103.0056, 140.1309}/*{117.4361,156.5599}*/;//M 为肤色在 YCbCr 颜色空间的样本均值 (Cb, Cr)，经验值 double C[2][2]={{160.1301,12.1430},//C 为肤色相似度模型的协方差矩阵，同上为经验值 {12.1430,299.4574}};// 注：因为运算仅需要该矩阵的逆矩阵值，故该值没有使用，仅作参考 double invC[2][2]={0.0077 ,-0.0041,-0.0041 ,0.0047 };//Ct 为 C的逆矩阵值，由 matlab计算而得 ////////////////////////////////////////////////////////////////////////// IplImage* pImg = img; double CrMean=0,CbMean=0,YMean=0; // 1 颜色转换：BGR->YCrCb IplImage*imgYCrCb=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U,3);// YCrCb图像 西南交通大学，信息学院，xgong@swjtu.cn 资料整理 5 cvCvtColor(pImg, imgYCrCb, CV_BGR2YCrCb);// 第 0,1,2层分别为 Y,Cr,Cb IplImage *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 根据 Sigle Gaussian Model计算颜色模型 IplImage *PCbCr=cvCreateImage(cvGetSize(pImg), IPL_DEPTH_32F, 1);//YCrCb颜色模型 IplImage *tempA=cvCreateImage(cvGetSize(pImg), IPL_DEPTH_32F, 1);//YCrCb颜色模型 IplImage *tempB=cvCreateImage(cvGetSize(pImg), IPL_DEPTH_32F, 1);//YCrCb颜色模型 cvSubS(imgCb32, cvScalar(M[0]), imgCb32);// x-m cvSubS(imgCr32, cvScalar(M[1]), imgCr32);// x-m cvAddWeighted(imgCb32, invC[0][0], imgCr32, invC[1][0], 0, tempA); cvAddWeighted(imgCb32, invC[0][1], imgCr32, invC[1][1], 0, tempB); cvMul(imgCb32, tempA, tempA, -0.5); cvMul(imgCr32, tempB, tempB, -0.5); cvAdd(tempA, tempB, PCbCr); cvExp(PCbCr, PCbCr); double max_val=0,min_val=0; cvMinMaxLoc(PCbCr,&min_val,&max_val); IplImage *proImg=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U, 1);//YCrCb颜色模型 double a=255/(max_val); cvConvertScaleAbs(PCbCr,proImg,a,0); m_proimg = cvCloneImage(proImg); if ((*skinImg)!=NULL) cvReleaseImage(skinImg); *skinImg = cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U, 1);//肤色结果 // 释放内存 cvReleaseImage(&proImg); cvReleaseImage(&imgYCrCb); cvReleaseImage(&imgY); cvReleaseImage(&imgCr); cvReleaseImage(&imgCb); cvReleaseImage(&imgY32); cvReleaseImage(&imgCr32); cvReleaseImage(&imgCb32); cvReleaseImage(&PCbCr); 西南交通大学，信息学院，xgong@swjtu.cn 资料整理 6 cvReleaseImage(&tempA); cvReleaseImage(&tempB); } 3.1高斯混合模型 （1）以下 matlab代码实现了高斯混合模型： function [Alpha, Mu, Sigma] = GMM_EM(Data, Alpha0, Mu0, Sigma0) %% EM 迭代停止条件 loglik_threshold = 1e-10; %% 初始化参数 [dim, N] = size(Data); M = size(Mu0,2); loglik_old = -realmax; nbStep = 0; Mu = Mu0; Sigma = Sigma0; Alpha = Alpha0; Epsilon = 0.0001; while (nbStep < 1200) nbStep = nbStep+1; %% E-步骤 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:M % PDF of each point Pxi(:,i) = GaussPDF(Data, Mu(:,i), Sigma(:,:,i)); end % 计算后验概率 beta(i|x) Pix_tmp = repmat(Alpha,[N 1]).*Pxi; Pix = Pix_tmp ./ (repmat(sum(Pix_tmp,2),[1 M])+realmin); Beta = sum(Pix); %% M-步骤 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:M % 更新权值 Alpha(i) = Beta(i) / N; % 更新均值 Mu(:,i) = Data*Pix(:,i) / Beta(i); % 更新方差 Data_tmp1 = Data - repmat(Mu(:,i),1,N); Sigma(:,:,i) = (repmat(Pix(:,i)',dim, 1) .* Data_tmp1*Data_tmp1') / Beta(i); %% Add a tiny variance to avoid numerical instability Sigma(:,:,i) = Sigma(:,:,i) + 1E-5.*diag(ones(dim,1)); end 西南交通大学，信息学院，xgong@swjtu.cn 资料整理 7 % %% Stopping criterion 1 %%%%%%%%%%%%%%%%%%%% % for i=1:M %Compute the new probability p(x|i) % Pxi(:,i) = GaussPDF(Data, Mu(:,i), Sigma(i)); % end %Compute the log likelihood % F = Pxi*Alpha'; % F(find(F