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

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

一、GMM 数学公式推导

1、预备知识：

（1）设离散型随机变量X 的分布率为： {} 2,1,P ===k p a X k k 则称()∑=

k

k k

p a

X E 为X 的数学期望或均值

（2）设连续型随机变量X 的概率密度函数（PDF ）为f(x) 其数学期望定义为：()()dx x xf X E ?

+∞

∞

-=

（3）()()()[]

2

X 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 n

T e

C

C x N 12

1

21

,;

其中，x 是维数为n 的样本向量（列向量），μ是期望，C 是协方差矩阵，|C|表示C 的行列式，1-C 表示C 的逆矩阵，()T

x μ-表示()μ-x 的转置。

3、混合高斯模型：GMM

设想有 m 个类：m 321????，，，， ，每类均服从正态分布。

各分布的中心点（均值）分别为：m 321μμμμ，，，，

方差分别为：m 321σσσσ，，，，

每一类在所有的类中所占的比例为 ()()()()m P P P P ????,,,,321 其中()11=∑=m

i i P ?。

同时，已知

个观察点：

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

示概率密度。

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

()()()()()()()()

()()()μμπ??σμ?σμ?σμ---

=-∑

=?++?+?=x C x m

i d i m m m T e

C

P P N P N P N x p 12

1

12221112,,,

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

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

可以直接采用MLE （最大似然估计）进行参数估计： 未知量为集合：()()()m P P C C ??μμλ,,1m 1m 1 ，，，，，，= 将衡量概率密度函数优劣的标准写出：()()∏==n

k k x P x p 1||λλ

即为：

()()

()

()()i k T i k x C x n k m

i d

i e

C P x p μμπ?λ---

==-∏∑

=12

1

11

|

|2|

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

大值的位置，就需用导数求极点，具体求解过程于下：

()()

()

()()

∑∑==∑==∏===n k P x N n

k x P x P x p m

i i i k k n

k k 1

,1|||1

1ln ln ln ln ?λλλλ

求导：()()()()()

()()()

()

()

()()}

|

|2{1

},{1ln ln 1-11

2

1

1

1

|1

1,1

,1|i k T i k k m

i i i k m

i i i k x C x m

i d

i m

k x p m

i i i k m

k P x N n

k P x N n

k x p e

C P P x N μμλ?λ?λλπ?λ?λλλλ---

======∑??=∑??∑=∑??

=??∑

∑

∑∑=

=

然后再分别对各个参数求导：①求参数i

μ ：②对 感兴趣，求偏导数有：③对

感兴趣，接下来的求导比较复杂，在此就没有继续推导。

（2）样本未知： EM 估计，算法流程： ①初始化：

方案1：协方差矩阵0j C 设为单位矩阵，每个模型比例的先验概率设为M

j 1

0=α，均值0j μ为随机数。

方案2：有K 均值（K-means)聚类算法对样本进行聚类，利用各类的均值作为

0j μ，并计算0j C ，0j α去各类样本占总数的比例。

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

()

M j n i x N x N M

k i k k

i j j ij ≤≤≤≤=

∑=1,1,||1

φα

φαβ

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

更新权值：n

n

i ij

j ∑==

1

β

α

更新均值：∑∑===

n i n

i ij

i j ij

x 1

1ββ

μ

更新方差矩阵：()()

∑∑==--=

n

i ij

n

i T

i

i

i

i

ij

j x x C 1

1

β

μμβ

④收敛条件：

不断地迭代步骤②和③，重复更新上面的三个值，直到()()εφφ<-||'|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 __,2

0?=σ ③

M

w 1

0=

④ 一般模型的个数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

σ

?i

最小的模型。新的高斯模型，用当前像素值作为新模型的均值，即()y x I i ,=μ，协方差为init std i _=σ，权重为

α=i w ，其中α为学习速率。

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

()()y x I t t i t i ,11,,αμαμ+-=- ⑥ ()()()2

1,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 #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, 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 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 + 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 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 + j*COMPONET + k] / 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])

{

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

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

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

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