混和高斯模型的推导和实现

1、基于GMM的运动目标检测方法研究1、 GMM数学公式推导1、 预备知识： （1）设离散型随机变量X的分布率为： 则称为X的数学期望或均值（2） 设连续型随机变量X的概率密度函数（PDF）为f(x) 其数学期望定义为：（3） 称为随机变量x的方差，称为X 的标准差（4） 正态分布： 概率密度函数为：（5） 设(x,y)为二维随机变量，若存在，则 称其为X和Y的协方差，记为cov(x,y) 2、 单高斯模型：SGM（也就是多维正态分布） 其概率密度函数PDF定义如下： 其中，x是维数为n的样本向量（列向量），是期望，C是协方差矩阵，|C|表示C的行列式，表示C的逆矩阵，表示的转置。3、 混合高斯模

2、型：GMM设想有 m个类： ，每类均服从正态分布。各分布的中心点（均值）分别为： 方差分别为：每一类在所有的类中所占的比例为 其中。同时，已知 个观察点： 。其中，用大写P表示概率，用小写p表示概率密度。则依此构想，可得概率密度函数为：其中d是维数，|·|是行列式但是在利用GMM进行目标检测时，这些模型的参数可能已知，也可能不知道，当参数已知时，可以直接利用GMM进行目标检测，在未知的情况下，需要对参数进行估计。对参数估计时，还要考虑样本分类是否已知。（1） 样本已知： 最大似然估计：可以直接采用MLE（最大似然估计）进行参数估计：未知量为集合： 将衡量概率密度函数优劣的标准写出：即

3、为：只要定出该标准的最大值位置，就可以求出最优的待定参数。为了 求出这个最大值的位置，就需用导数求极点，具体求解过程于下：求导：然后再分别对各个参数求导：求参数 ：对 感兴趣，求偏导数有：对 感兴趣，接下来的求导比较复杂，在此就没有继续推导。（2） 样本未知： EM估计，算法流程：初始化： 方案1：协方差矩阵设为单位矩阵，每个模型比例的先验概率设为，均值为随机数。 方案2：有K均值（K-means)聚类算法对样本进行聚类，利用各类的均值作为，并计算，去各类样本占总数的比例。估计步骤（E-step）：令的后验概率为： 最大化步骤（M-step）：更新权值：更新均值：更新方差矩阵：收敛条件： 不断

4、地迭代步骤和，重复更新上面的三个值，直到，其中为更新参数后计算的值，即前后两次迭代得到的结果变化小于一定程度则终止迭代，通常2、 GMM发展历史及现状背景建模方法有很多种，如中值法、均值法、卡尔曼滤波器模型、码本背景模型等，其中混合高斯模型是最经典的算法。GMM最早是由CHris Stauffer等在1中提出的，该方法是按照高斯分布对每个像素建立模型， 并通过基于回归滤波的在线 EM 近似方法对模型参数进行更新，它能鲁棒地克服光照变化、 树枝摇动等造成的影响，但该方法也存在一些问题：1）该方法对运动物体在场景中停止不动或者长时间停止时检测失效，而且带有初始学习速度慢，在线更新费时、计算量大；2

5、）无法完整准确地检测大并且运动缓慢的运动目标，运动目标的像素点不集中，只能检测到运动目标的部分轮廓，无法提取出目标对象的完整区域；3）无法将背景显露区域与运动目标区域很好地区分开；4）当运动目标由静止缓慢转化为运动时，易将背景显露区检测为前景，出现“影子”现象。3、 GMM缺点及改进方法针对上述问题，一些科学研究者又在GMM算法的基础上做了很多的改进：张、白等人2引入分块思想，把图像分为L*L块；黄、胡等人3也引入了分块的思想，但是他们的分块理念是以当前像素点的8邻域作为一块；华、刘4把GMM与改进的帧差法(相邻两帧图像对应像素点8邻域像素值相减之和)相结合，提高了计算效率；Suo等人5是将混

6、合高斯模型中的模型个数采改进为自适应的；刘等人6融合帧间差分法，检测背景显露区域和运动区域，很好的解决了问题4。除此之外，还有基于纹理的混合高斯模型。4、 GMM算法流程（1） 用第一帧图像对高斯混合模型进行初始化 一般模型的个数M为3-6个，其中std_init设置为20（2） 对于t时刻的像素，分别与已经存在的M个高斯模型依次进行匹配： （3） 如果满足匹配条件，则该像素值与高斯模型匹配成功。如果匹配不成功： a：当k<K时，增加新的高斯模型； b：当k=K时，用新高斯模型代替优先级最小的模型。新的高斯模型，用当前像素值作为新模型的均值，即，协方差为，权重为，其中为学习速率。（4）

7、未匹配模式的均值和方差不变，对匹配模式的第i个高斯模型参数进行更新： （5） 高斯模型参数更新完毕后，对每个像素点的K歌高斯模型按优先级降序排序。取前B个高斯模型作为背景像素的最佳描述： （6） 继续对与上述B个高斯模型进行匹配检验，如果与前B个高斯模型的任意一个匹配，则该像素点为背景点；否则为前景点。（7） 重复步骤（2）-（6），直到视频结束。5、 GMM代码实现#include<opencv.hpp>#include<highgui.h>#include<cv.h>using namespace cv;using namespace std;#defi

8、ne COMPONET 5 /混合高斯模型个数#define ALPHA 0.03 /学习率#define SD_INIT 6 /方差初始值#define THRESHOLD 0.25 /前景所占比例#define D 2.5 int main()CvCapture *capture = cvCreateFileCapture("E:project2videosvideo.avi");IplImage *frame, *grayFrame, *foreground, *background;int *foreg, *backg, *rank_index;double *we

9、ight, *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 = cv

10、CreateImage(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)

11、*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 means sigma = (double *)malloc(sizeof(double)*width*height*COMPONET); /pixel standard dev

12、iations u_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+)meani*width*COMPONET + j*COMPONET + k = cvRandReal(&state) * 255;sigmai*

13、width*COMPONET + j*COMPONET + k = SD_INIT;weighti*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 mean for (int i = 0; i < height; i+)for (int j =

14、 0; j < width; j+)for (int k = 0; k < COMPONET; k+)u_diffi*width*COMPONET + j*COMPONET + k = abs(uchar)grayFrame->imageDatai*width + j - meani*width*COMPONET + j*COMPONET + k);/update gaussian components for each pixel for (int i = 0; i < height; i+)for (int j = 0; j < width; j+)match

15、 = 0;double sum_weight = 0;for (int k = 0; k < COMPONET; k+)if (u_diffi*width*COMPONET + j*COMPONET + k <= D*sigmai*width*COMPONET + j*COMPONET + k) /pixel matches component match = 1; / variable to signal component match /update weights, mean, sd, p weighti*width*COMPONET + j*COMPONET + k = (

16、1 - ALPHA)*weighti*width*COMPONET + j*COMPONET + k + ALPHA;/*p = ALPHA / weighti*width*COMPONET + j*COMPONET + k;meani*width*COMPONET + j*COMPONET + k = (1 - p)*meani*width*COMPONET + j*COMPONET + k + p*(uchar)grayFrame->imageDatai*width + j;sigmai*width*COMPONET + j*COMPONET + k = sqrt(1 - p)*(s

17、igmai*width*COMPONET + j*COMPONET + k * sigmai*width*COMPONET + j*COMPONET + k) + p*(pow(uchar)grayFrame->imageDatai*width + j - meani*width*COMPONET + j*COMPONET + k, 2);*/meani*width*COMPONET + j*COMPONET + k = (1 - ALPHA)*meani*width*COMPONET + j*COMPONET + k + ALPHA*(uchar)grayFrame->image

18、Datai*width + j;sigmai*width*COMPONET + j*COMPONET + k = sqrt(1 - ALPHA)*(sigmai*width*COMPONET + j*COMPONET + k * sigmai*width*COMPONET + j*COMPONET + k) + ALPHA*(pow(uchar)grayFrame->imageDatai*width + j - meani*width*COMPONET + j*COMPONET + k, 2);/else/weighti*width*COMPONET + j*COMPONET + k =

19、 (1 - ALPHA)*weighti*width*COMPONET + j*COMPONET + k; / weight slighly decreases /sum_weight += weighti*width*COMPONET + j*COMPONET + k;/权重归一化for (int k = 0; k < COMPONET; k+)weighti*width*COMPONET + j*COMPONET + k = weighti*width*COMPONET + j*COMPONET + k / sum_weight;/获取权重最小下标double temp = weig

20、hti*width*COMPONET + j*COMPONET;int min_index = 0;backgi*width + j = 0;for (int k = 0; k < COMPONET; k+)backgi*width + j = backgi*width + j + meani*width*COMPONET + j*COMPONET + k * weighti*width*COMPONET + j*COMPONET + k;if (weighti*width*COMPONET + j*COMPONET + k < temp)min_index = k;temp =

21、weighti*width*COMPONET + j*COMPONET + k;rank_indexk = k;background->imageDatai*width + j = (uchar)backgi*width + j;/if no components match, create new component if (match = 0)meani*width*COMPONET + j*COMPONET + min_index = (uchar)grayFrame->imageDatai*width + j;sigmai*width*COMPONET + j*COMPON

22、ET + min_index = SD_INIT;weighti*width*COMPONET + j*COMPONET + min_index = 1 / COMPONET;/计算优先级for (int k = 0; k < COMPONET; k+)rankk = weighti*width*COMPONET + j*COMPONET + k / sigmai*width*COMPONET + j*COMPONET + k;/sort rank values for (int k = 1; k < COMPONET; k+)for (int m = 0; m < k; m

23、+)if (rankk > rankm)/swap max values rank_temp = rankm;rankm = rankk;rankk = rank_temp;/swap max index values rank_index_temp = rank_indexm;rank_indexm = rank_indexk;rank_indexk = rank_index_temp;/calculate foreground match = 0;int b = 0; while (match = 0) && (b < COMPONET)if (weighti*width*COMPONET + j*COMPONET + rank_indexb >= THRESHOLD