基于GMM 的运动目标检测方法研究一、GMM 数学公式推导1、预备知识：（1）设离散型随机变量X 的分布率为： {} 2,1,P ===k p a X k k 则称()∑=kk kp aX E 为X 的数学期望或均值（2）设连续型随机变量X 的概率密度函数（PDF ）为f(x) 其数学期望定义为：()()dx x xf X E ⎰+∞∞-=（3）()()()[]2X E X E X D -=称为随机变量x 的方差，()X D 称为X的标准差（4）正态分布：()2,~σμN X 概率密度函数为：()()⎥⎥⎦⎤⎢⎢⎣⎡--=22221σμσπx e x p（5）设(x,y)为二维随机变量，()[]()[]{}Y E Y X E X E --若存在，则 称其为X 和Y 的协方差，记为cov(x,y)()()[]()[]{}()XY E Y E Y X E X E Y X =--=,cov 2、单高斯模型：SGM （也就是多维正态分布） 其概率密度函数PDF 定义如下： ()()()()μμπμ----=x C x nT eCC x N 12121,;其中，x 是维数为n 的样本向量（列向量），μ是期望，C 是协方差矩阵，|C|表示C 的行列式，1-C 表示C 的逆矩阵，()Tx μ-表示()μ-x 的转置。

3、混合高斯模型：GMM设想有 m 个类：m 321ϖϖϖϖ，，，， ，每类均服从正态分布。

各分布的中心点（均值）分别为：m 321μμμμ，，，，方差分别为：m 321σσσσ，，，，每一类在所有的类中所占的比例为 ()()()()m P P P P ϖϖϖϖ,,,,321 其中()11=∑=mi i P ϖ。

同时，已知个观察点：。

其中，用大写P 表示概率，用小写p 表示概率密度。

则依此构想，可得概率密度函数为：()()()()()()()()()()()μμπϖϖσμϖσμϖσμ---=-∑=⋅++⋅+⋅=x C x mi d i m m m T eCP P N P N P N x p 12112221112,,,其中d 是维数，|·|是行列式但是在利用GMM 进行目标检测时，这些模型的参数可能已知，也可能不知道，当参数已知时，可以直接利用GMM 进行目标检测，在未知的情况下，需要对参数进行估计。

对参数估计时，还要考虑样本分类是否已知。

（1）样本已知： 最大似然估计：可以直接采用MLE （最大似然估计）进行参数估计： 未知量为集合：()()()m P P C C ϖϖμμλ,,1m 1m 1 ，，，，，，= 将衡量概率密度函数优劣的标准写出：()()∏==nk k x P x p 1||λλ即为：()()()()()i k T i k x C x n k mi di eC P x p μμπϖλ---==-∏∑=12111||2|只要定出该标准的最大值位置，就可以求出最优的待定参数。

为了 求出这个最大值的位置，就需用导数求极点，具体求解过程于下：()()()()()∑∑==∑==∏===n k P x N nk x P x P x p mi i i k k nk k 1,1|||11ln ln ln ln ϖλλλλ求导：()()()()()()()()()()()()}||2{1},{1ln ln 1-112111|11,1,1|i k T i k k mi i i k mi i i k x C x mi di mk x p mi i i k mk P x N nk P x N nk x p eC P P x N μμλϖλϖλλπϖλϖλλλλ---======∑∂∂=∑∂∂∑=∑∂∂=∂∂∑∑∑∑==然后再分别对各个参数求导：①求参数iμ ：②对 感兴趣，求偏导数有：③对感兴趣，接下来的求导比较复杂，在此就没有继续推导。

（2）样本未知： EM 估计，算法流程： ①初始化：方案1：协方差矩阵0j C 设为单位矩阵，每个模型比例的先验概率设为Mj 10=α，均值0j μ为随机数。

方案2：有K 均值（K-means)聚类算法对样本进行聚类，利用各类的均值作为0j μ，并计算0j C ，0j α去各类样本占总数的比例。

②估计步骤（E-step ）： 令j α的后验概率为： ()()M j n i x N x N Mk i k ki j j ij ≤≤≤≤=∑=1,1,||1φαφαβ③最大化步骤（M-step ）：更新权值：nni ijj ∑==1βα更新均值：∑∑===n i ni iji j ijx 11ββμ更新方差矩阵：()()∑∑==--=ni ijni Tiiiiijj x x C 11βμμβ④收敛条件：不断地迭代步骤②和③，重复更新上面的三个值，直到()()εφφ<-||'|X p X p ，其中为更新参数后计算的值，即前后两次迭代得到的结果变化小于一定程度则终止迭代，通常-510=ε 二、GMM 发展历史及现状背景建模方法有很多种，如中值法、均值法、卡尔曼滤波器模型、码本背景模型等，其中混合高斯模型是最经典的算法。

GMM 最早是由CHris Stauffer 等在[1]中提出的，该方法是按照高斯分布对每个像素建立模型， 并通过基于回归滤波的在线 EM 近似方法对模型参数进行更新，它能鲁棒地克服光照变化、 树枝摇动等造成的影响，但该方法也存在一些问题：1）该方法对运动物体在场景中停止不动或者长时间停止时检测失效，而且带有初始学习速度慢，在线更新费时、计算量大；2）无法完整准确地检测大并且运动缓慢的运动目标，运动目标的像素点不集中，只能检测到运动目标的部分轮廓，无法提取出目标对象的完整区域；3）无法将背景显露区域与运动目标区域很好地区分开；4）当运动目标由静止缓慢转化为运动时，易将背景显露区检测为前景，出现“影子”现象。

三、GMM 缺点及改进方法针对上述问题，一些科学研究者又在GMM 算法的基础上做了很多的改进：张、白等人[2]引入分块思想，把图像分为L*L 块；黄、胡等人[3]也引入了分块的思想，但是他们的分块理念是以当前像素点的8邻域作为一块；华、刘[4]把GMM 与改进的帧差法(相邻两帧图像对应像素点8邻域像素值相减之和)相结合，提高了计算效率；Suo 等人[5]是将混合高斯模型中的模型个数采改进为自适应的；刘等人[6]融合帧间差分法，检测背景显露区域和运动区域，很好的解决了问题4。

除此之外，还有基于纹理的混合高斯模型。

四、GMM 算法流程（1）用第一帧图像对高斯混合模型进行初始化 ()()0,,,0y x I y x =μ ① ()init std y x _,0=σ ②()init std init std y x __,20⨯=σ ③Mw 10=④ 一般模型的个数M 为3-6个，其中std_init 设置为20（2）对于t 时刻的像素()y x I t ,，分别与已经存在的M 个高斯模型依次进行匹配：()1,1,5.2|),(,--<-t i t i t y x y x I σμ ⑤（3）如果满足匹配条件，则该像素值与高斯模型匹配成功。

如果匹配不成功： a ：当k<K 时，增加新的高斯模型； b ：当k=K 时，用新高斯模型代替优先级σϖi最小的模型。

新的高斯模型，用当前像素值作为新模型的均值，即()y x I i ,=μ，协方差为init std i _=σ，权重为α=i w ，其中α为学习速率。

（4）未匹配模式的均值和方差不变，对匹配模式的第i 个高斯模型参数进行更新：()()y x I t t i t i ,11,,αμαμ+-=- ⑥ ()()()21,21,2,,1---+-=t i t t i t i y x I μασασ ⑦()αα+-=-1,,1t i t i w w ⑧（5）高斯模型参数更新完毕后，对每个像素点的K 歌高斯模型按优先级σϖi降序排序。

取前B 个高斯模型作为背景像素的最佳描述：15.0;min arg 1<<⎪⎭⎫⎝⎛>=∑=T T w B M k i k⑨（6）继续对()y x I t ,与上述B 个高斯模型进行匹配检验，如果()y x I t ,与前B 个高斯模型的任意一个匹配，则该像素点为背景点；否则为前景点。

（7）重复步骤（2）-（6），直到视频结束。

五、GMM 代码实现 #include<opencv.hpp> #include<highgui.h> #include<cv.h>using namespace cv; using namespace std;#define COMPONET 5 //混合高斯模型个数 #define ALPHA 0.03 //学习率 #define SD_INIT 6 //方差初始值 #define THRESHOLD 0.25 //前景所占比例 #define D 2.5int main() {CvCapture*capture=cvCreateFileCapture("E:\\project2\\videos\\video.avi");IplImage *frame, *grayFrame, *foreground, *background;int *foreg, *backg, *rank_index;double *weight, *mean, *sigma, *u_diff, *rank;double p = ALPHA / (1 / (double)COMPONET);double rank_temp = 0;int rank_index_temp = 0;CvRNG state; //随机生成状态器int match, height, width;frame = cvQueryFrame(capture);grayFrame = cvCreateImage(CvSize(frame->width, frame->height), IPL_DEPTH_8U, 1);foreground = cvCreateImage(CvSize(frame->width, frame->height), IPL_DEPTH_8U, 1);background = cvCreateImage(CvSize(frame->width, frame->height), IPL_DEPTH_8U, 1);height = grayFrame->height;width = grayFrame->widthStep;foreg = (int*)malloc(sizeof(int)*width*height);backg = (int*)malloc(sizeof(int)*width*height);rank = (double*)malloc(sizeof(double) * 1 * COMPONET); //优先级weight = (double*)malloc(sizeof(double)*width*height*COMPONET); //权重mean = (double *)malloc(sizeof(double)*width*height*COMPONET);//pixel meanssigma = (double *)malloc(sizeof(double)*width*height*COMPONET);//pixel standard deviationsu_diff = (double *)malloc(sizeof(double)*width*height*COMPONET);//difference of each pixel from mean//初始化均值、方差、权重for (int i = 0; i < height; i++){for (int j = 0; j < width; j++){for (int k = 0; k < COMPONET; k++){mean[i*width*COMPONET + j*COMPONET + k] = cvRandReal(&state) * 255;sigma[i*width*COMPONET + j*COMPONET + k] = SD_INIT;weight[i*width*COMPONET + j*COMPONET + k] = (double)1 / COMPONET;}}}while (1){rank_index = (int *)malloc(sizeof(int)*COMPONET);cvCvtColor(frame, grayFrame, CV_BGR2GRAY);// calculate difference of pixel values from meanfor (int i = 0; i < height; i++){for (int j = 0; j < width; j++){for (int k = 0; k < COMPONET; k++){u_diff[i*width*COMPONET + j*COMPONET + k] =abs((uchar)grayFrame->imageData[i*width + j] -mean[i*width*COMPONET + j*COMPONET + k]);}}}//update gaussian components for each pixelfor (int i = 0; i < height; i++){for (int j = 0; j < width; j++){match = 0;double sum_weight = 0;for (int k = 0; k < COMPONET; k++){if (u_diff[i*width*COMPONET + j*COMPONET + k]<= D*sigma[i*width*COMPONET + j*COMPONET + k]) //pixelmatches component{match = 1;// variable to signal component match//update weights, mean, sd, pweight[i*width*COMPONET + j*COMPONET + k] = (1 - ALPHA)*weight[i*width*COMPONET + j*COMPONET + k] + ALPHA;/*p = ALPHA / weight[i*width*COMPONET + j*COMPONET + k];mean[i*width*COMPONET + j*COMPONET + k] = (1 - p)*mean[i*width*COMPONET + j*COMPONET + k] + p*(uchar)grayFrame->imageData[i*width + j];sigma[i*width*COMPONET + j*COMPONET + k] = sqrt((1 - p)*(sigma[i*width*COMPONET + j*COMPONET + k] * sigma[i*width*COMPONET + j*COMPONET + k]) + p*(pow((uchar)grayFrame->imageData[i*width + j] - mean[i*width*COMPONET + j*COMPONET + k], 2)));*/mean[i*width*COMPONET + j*COMPONET + k] = (1 - ALPHA)*mean[i*width*COMPONET + j*COMPONET + k] + ALPHA*(uchar)grayFrame->imageData[i*width + j];sigma[i*width*COMPONET + j*COMPONET + k] = sqrt((1 - ALPHA)*(sigma[i*width*COMPONET + j*COMPONET + k] * sigma[i*width*COMPONET + j*COMPONET + k]) + ALPHA*(pow((uchar)grayFrame->imageData[i*width + j] -mean[i*width*COMPONET + j*COMPONET + k], 2)));}//else{// weight[i*width*COMPONET + j*COMPONET + k] =(1 - ALPHA)*weight[i*width*COMPONET + j*COMPONET + k]; // weight slighly decreases//}sum_weight += weight[i*width*COMPONET +j*COMPONET + k];}//权重归一化for (int k = 0; k < COMPONET; k++){weight[i*width*COMPONET + j*COMPONET + k] =weight[i*width*COMPONET + j*COMPONET + k] / sum_weight;}//获取权重最小下标double temp = weight[i*width*COMPONET +j*COMPONET];int min_index = 0;backg[i*width + j] = 0;for (int k = 0; k < COMPONET; k++){backg[i*width + j] = backg[i*width + j] + mean[i*width*COMPONET + j*COMPONET + k] * weight[i*width*COMPONET + j*COMPONET + k];if (weight[i*width*COMPONET + j*COMPONET + k] < temp){min_index = k;temp = weight[i*width*COMPONET + j*COMPONET + k];}rank_index[k] = k;}background->imageData[i*width + j] = (uchar)backg[i*width + j];//if no components match, create new componentif (match == 0){mean[i*width*COMPONET + j*COMPONET + min_index]= (uchar)grayFrame->imageData[i*width + j];sigma[i*width*COMPONET + j*COMPONET + min_index] = SD_INIT;weight[i*width*COMPONET + j*COMPONET + min_index] = 1 / COMPONET;}//计算优先级for (int k = 0; k < COMPONET; k++){rank[k] = weight[i*width*COMPONET + j*COMPONET + k] / sigma[i*width*COMPONET + j*COMPONET + k];}//sort rank valuesfor (int k = 1; k < COMPONET; k++){for (int m = 0; m < k; m++){if (rank[k] > rank[m]){//swap max valuesrank_temp = rank[m];rank[m] = rank[k];rank[k] = rank_temp;//swap max index valuesrank_index_temp = rank_index[m];rank_index[m] = rank_index[k];rank_index[k] = rank_index_temp;}}}//calculate foregroundmatch = 0;int b = 0;while ((match == 0) && (b < COMPONET)){if (weight[i*width*COMPONET + j*COMPONET + rank_index[b]] >= THRESHOLD){if (abs(u_diff[i*width*COMPONET + j*COMPONET + rank_index[b]]) <= D*sigma[i*width*COMPONET + j*COMPONET + rank_index[b]]){foreground->imageData[i*width + j] = 0;match = 1;}else{foreground->imageData[i*width + j] = (uchar)grayFrame->imageData[i*width + j];}}b++;}}}frame = cvQueryFrame(capture);cvShowImage("fore", foreground);cvShowImage("back", background);cvShowImage("frame", frame);char s = cvWaitKey(33);if (s == 27) break;free(rank_index);}return 0;}六、参考文献[1]Chris Stauffer,W.E.L Grimson.Adaptive background mixture models for real-time tracking[2]张燕平、白云球.应用改进混合高斯模型的运动目标检测[3]黄大卫、胡文翔。