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

基于GMM的运动目标检测方法研究

一、GMM数学公式推导

1、预备知识：

（1）设离散型随机变量X的分布率为：

 PXakpk,k1,2

EXakpkk则称为X的数学期望或均值

（2）设连续型随机变量X的概率密度函数（PDF）为f(x)

其数学期望定义为：

EXxfxdx

2DXEXEX（3）称为随机变量x的方差，DX称为X

的标准差

2 X~N,（4）正态分布：

概率密度函数为：

px1e2x222

（5）设(x,y)为二维随机变量，EXEXYEY若存在，则

称其为X和Y的协方差，记为cov(x,y)

covX,YEXEXYEYEXY 2、单高斯模型：SGM（也就是正态分布）

其概率密度函数PDF定义如下：

Nx;,C

12nCe1xTC1x2

其中，x是维数为n的样本向量（列向量），μ是期望，C是协方差矩阵，|C|表示C

1的行列式，C表示C的逆矩阵，x表示x的转置。

T3、混合高斯模型：GMM

，m设想有 m个类：1，2，3， ，每类均服从正态分布。

，m各分布的中心点（均值）分别为：1，2，3，

，m 方差分别为：1，2，3，

每一类在所有的类中所占的比例为 P1,P2,P3,,Pm

其中

P1ii1m。

同时，已知 个观察点： 密度。

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

则依此构想，可得概率密度函数为：

pxN1,1P1N2,2P2Nm,mPmi1mPi2dCe1xTC1x2

其中d是维数，|·|是行列式

但是在利用GMM进行目标检测时，这些模型的参数可能已知，也可能不知道，当参数已知时，可以直接利用GMM进行目标检测，在未知的情况下，需要对参数进行估计。对参数估计时，还要考虑样本分类是否已知。

（1）样本已知：

最大似然估计：

可以直接采用MLE（最大似然估计）进行参数估计：

，m，C1，,Cm，P1，,Pm未知量为集合：1，

将衡量概率密度函数优劣的标准写出：

px|Pxk|k1n

即为：

px|k1i1nmPi2d|C|e1xkiTC1xki2

只要定出该标准的最大值位置，就可以求出最优的待定参数。为了 求出这个最大值的

lnpx|lnk1lnPxk|k1nnPxk|n位置，就需用导数求极点，具体求解过程于下：

lni1k1Nxk,iPi

mNxk,iPipx|lnlni1k1k1nnmk1mm1Nxk,iPii1mm{Nxk,iPi}i1Pie1xkiTC-1xki2求导：

k11pxk|m{i12d|C|}

然后再分别对各个参数求导：①求参数i ：②对 感兴趣，求偏导数有：③对 感兴趣，接下来的求导比较复杂，在此就没有继续推导。

（2）样本未知：

EM估计，算法流程：

①初始化：

方案1：协方差矩阵Cj0设为单位矩阵，每个模型比例的先验概率设为为随机数。

j01M，均值j0

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

②估计步骤（E-step）： 令j的后验概率为：

jNjxi|ij

k1MkNkxi|,1in,1jM

③最大化步骤（M-step）：

更新权值：

ji1nijn

j

xii1

ni1

n

ij

更新均值：

ij

n

Cjxxijiiiii1T更新方差矩阵：

i1nij

④收敛条件：

不断地迭代步骤②和③，重复更新上面的三个值，直到pX|p'X||，其中为更

-5新参数后计算的值，即前后两次迭代得到的结果变化小于一定程度则终止迭代，通常10

二、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x,yIx,y,0 ① 0x,ystd_init ②

02x,ystd_initstd_init ③

w01M ④

一般模型的个数M为3-6个，其中std_init设置为20

（2）对于t时刻的像素Itx,y，分别与已经存在的M个高斯模型依次进行匹配：

Itx,yi,t1(x,y)|2.5i,t1 ⑤

（3）如果满足匹配条件，则该像素值与高斯模型匹配成功。如果匹配不成功：

a：当kib：当k=K时，用新高斯模型代替优先级最小的模型。新的高斯模型，用当前像素

值作为新模型的均值，即iIx,y，协方差为istd_init，权重为wi，其中α为学习速率。

（4）未匹配模式的均值和方差不变，对匹配模式的第i个高斯模型参数进行更新：

i,t1i,t1Itx,y ⑥

i2,t1i2,t1Itx,yi,t12 ⑦

wi,t1wi,t1 ⑧

i（5）高斯模型参数更新完毕后，对每个像素点的K歌高斯模型按优先级降序排序。

取前B个高斯模型作为背景像素的最佳描述：

MBargminwiT;0.5T1kk1 ⑨

（6）继续对Itx,y与上述B个高斯模型进行匹配检验，如果Itx,y与前B个高斯模型的任意一个匹配，则该像素点为背景点；否则为前景点。

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

五、GMM代码实现

#include

#include

#include

using namespace cv;

using namespace std;

#define COMPONET 5 //混合高斯模型个数

#define ALPHA 0.03 //学习率

#define SD_INIT 6 //方差初始值

#define THRESHOLD 0.25 //前景所占比例

#define D 2.5

int 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, IPL_DEPTH_8U, 1);

foreground = cvCreateImage(CvSize(frame->width, IPL_DEPTH_8U, 1);

background = cvCreateImage(CvSize(frame->width, IPL_DEPTH_8U, 1);

height = grayFrame->height;

width = grayFrame->widthStep;

foreg = (int*)malloc(sizeof(int)*width*height);

backg = (int*)malloc(sizeof(int)*width*height);

frame->height),

frame->height),

frame->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 deviations

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++)

{

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 mean

for (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 pixel for (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]) //pixel matches component

{

match = 1; // variable to signal component match

//update weights, mean, sd, p

weight[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 p)*mean[i*width*COMPONET

+ +

j*COMPONET + k] +

= k]

(1 - +

j*COMPONET

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 k]

- +

ALPHA)*mean[i*width*COMPONET j*COMPONET

ALPHA*(uchar)grayFrame->imageData[i*width + j];

sigma[i*width*COMPONET + j*COMPONET +

+ k] = +

sqrt((1 k] k]) j]

- * + -

ALPHA)*(sigma[i*width*COMPONET sigma[i*width*COMPONET

+

j*COMPONET

+ +

j*COMPONET

ALPHA*(pow((uchar)grayFrame->imageData[i*width 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 component

if (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 sigma[i*width*COMPONET + j*COMPONET + k];

}

//sort rank values

for (int k = 1; k < COMPONET; k++)

{

for (int m = 0; m < k; m++)

{

if (rank[k] > rank[m])

{

j*COMPONET k] /

+ +

//swap max values

rank_temp = rank[m];

rank[m] = rank[k];

rank[k] = rank_temp;

//swap max index values

rank_index_temp = rank_index[m];

rank_index[m] = rank_index[k];

rank_index[k] = rank_index_temp;

}

}

}

//calculate foreground

match = 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\

cvShowImage(\"back\

cvShowImage(\"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]黄大卫、胡文翔。改进单高斯模型的视频前景提取与破碎目标合并算法

[4]华媛蕾、刘万军.改进混合高斯模型的运动目标检测算法

[5]Peng Suo, Yanjiang Wang.Improved Adaptive Gaussian Mixture Model for Background Subtraction

[6]刘鑫、刘辉.混合高斯模型和帧间差分相融合的自适应背景模型