混和高斯模型的推導(dǎo)和實(shí)現(xiàn)

1、基于GMM的運(yùn)動(dòng)目標(biāo)檢測(cè)方法研究1、 GMM數(shù)學(xué)公式推導(dǎo)1、 預(yù)備知識(shí)： （1）設(shè)離散型隨機(jī)變量X的分布率為： 則稱為X的數(shù)學(xué)期望或均值（2） 設(shè)連續(xù)型隨機(jī)變量X的概率密度函數(shù)（PDF）為f(x) 其數(shù)學(xué)期望定義為：（3） 稱為隨機(jī)變量x的方差，稱為X 的標(biāo)準(zhǔn)差（4） 正態(tài)分布： 概率密度函數(shù)為：（5） 設(shè)(x,y)為二維隨機(jī)變量，若存在，則 稱其為X和Y的協(xié)方差，記為cov(x,y) 2、 單高斯模型：SGM（也就是多維正態(tài)分布） 其概率密度函數(shù)PDF定義如下： 其中，x是維數(shù)為n的樣本向量（列向量），是期望，C是協(xié)方差矩陣，|C|表示C的行列式，表示C的逆矩陣，表示的轉(zhuǎn)置。3、 混合高斯模

2、型：GMM設(shè)想有 m個(gè)類： ，每類均服從正態(tài)分布。各分布的中心點(diǎn)（均值）分別為： 方差分別為：每一類在所有的類中所占的比例為 其中。同時(shí)，已知 個(gè)觀察點(diǎn)： 。其中，用大寫P表示概率，用小寫p表示概率密度。則依此構(gòu)想，可得概率密度函數(shù)為：其中d是維數(shù)，|是行列式但是在利用GMM進(jìn)行目標(biāo)檢測(cè)時(shí)，這些模型的參數(shù)可能已知，也可能不知道，當(dāng)參數(shù)已知時(shí)，可以直接利用GMM進(jìn)行目標(biāo)檢測(cè)，在未知的情況下，需要對(duì)參數(shù)進(jìn)行估計(jì)。對(duì)參數(shù)估計(jì)時(shí)，還要考慮樣本分類是否已知。（1） 樣本已知： 最大似然估計(jì)：可以直接采用MLE（最大似然估計(jì)）進(jìn)行參數(shù)估計(jì)：未知量為集合： 將衡量概率密度函數(shù)優(yōu)劣的標(biāo)準(zhǔn)寫出：即為：只要定出該

3、標(biāo)準(zhǔn)的最大值位置，就可以求出最優(yōu)的待定參數(shù)。為了 求出這個(gè)最大值的位置，就需用導(dǎo)數(shù)求極點(diǎn)，具體求解過程于下：求導(dǎo)：然后再分別對(duì)各個(gè)參數(shù)求導(dǎo)：求參數(shù) ：對(duì) 感興趣，求偏導(dǎo)數(shù)有：對(duì) 感興趣，接下來的求導(dǎo)比較復(fù)雜，在此就沒有繼續(xù)推導(dǎo)。（2） 樣本未知： EM估計(jì)，算法流程：初始化： 方案1：協(xié)方差矩陣設(shè)為單位矩陣，每個(gè)模型比例的先驗(yàn)概率設(shè)為，均值為隨機(jī)數(shù)。 方案2：有K均值（K-means)聚類算法對(duì)樣本進(jìn)行聚類，利用各類的均值作為，并計(jì)算，去各類樣本占總數(shù)的比例。估計(jì)步驟（E-step）：令的后驗(yàn)概率為： 最大化步驟（M-step）：更新權(quán)值：更新均值：更新方差矩陣：收斂條件： 不斷地迭代步驟和，

4、重復(fù)更新上面的三個(gè)值，直到，其中為更新參數(shù)后計(jì)算的值，即前后兩次迭代得到的結(jié)果變化小于一定程度則終止迭代，通常2、 GMM發(fā)展歷史及現(xiàn)狀背景建模方法有很多種，如中值法、均值法、卡爾曼濾波器模型、碼本背景模型等，其中混合高斯模型是最經(jīng)典的算法。GMM最早是由CHris Stauffer等在1中提出的，該方法是按照高斯分布對(duì)每個(gè)像素建立模型， 并通過基于回歸濾波的在線 EM 近似方法對(duì)模型參數(shù)進(jìn)行更新，它能魯棒地克服光照變化、 樹枝搖動(dòng)等造成的影響，但該方法也存在一些問題：1）該方法對(duì)運(yùn)動(dòng)物體在場(chǎng)景中停止不動(dòng)或者長(zhǎng)時(shí)間停止時(shí)檢測(cè)失效，而且?guī)в谐跏紝W(xué)習(xí)速度慢，在線更新費(fèi)時(shí)、計(jì)算量大；2）無法完整準(zhǔn)確

5、地檢測(cè)大并且運(yùn)動(dòng)緩慢的運(yùn)動(dòng)目標(biāo)，運(yùn)動(dòng)目標(biāo)的像素點(diǎn)不集中，只能檢測(cè)到運(yùn)動(dòng)目標(biāo)的部分輪廓，無法提取出目標(biāo)對(duì)象的完整區(qū)域；3）無法將背景顯露區(qū)域與運(yùn)動(dòng)目標(biāo)區(qū)域很好地區(qū)分開；4）當(dāng)運(yùn)動(dòng)目標(biāo)由靜止緩慢轉(zhuǎn)化為運(yùn)動(dòng)時(shí)，易將背景顯露區(qū)檢測(cè)為前景，出現(xiàn)“影子”現(xiàn)象。3、 GMM缺點(diǎn)及改進(jìn)方法針對(duì)上述問題，一些科學(xué)研究者又在GMM算法的基礎(chǔ)上做了很多的改進(jìn)：張、白等人2引入分塊思想，把圖像分為L(zhǎng)*L塊；黃、胡等人3也引入了分塊的思想，但是他們的分塊理念是以當(dāng)前像素點(diǎn)的8鄰域作為一塊；華、劉4把GMM與改進(jìn)的幀差法(相鄰兩幀圖像對(duì)應(yīng)像素點(diǎn)8鄰域像素值相減之和)相結(jié)合，提高了計(jì)算效率；Suo等人5是將混合高斯模型中的

6、模型個(gè)數(shù)采改進(jìn)為自適應(yīng)的；劉等人6融合幀間差分法，檢測(cè)背景顯露區(qū)域和運(yùn)動(dòng)區(qū)域，很好的解決了問題4。除此之外，還有基于紋理的混合高斯模型。4、 GMM算法流程（1） 用第一幀圖像對(duì)高斯混合模型進(jìn)行初始化 一般模型的個(gè)數(shù)M為3-6個(gè)，其中std_init設(shè)置為20（2） 對(duì)于t時(shí)刻的像素，分別與已經(jīng)存在的M個(gè)高斯模型依次進(jìn)行匹配： （3） 如果滿足匹配條件，則該像素值與高斯模型匹配成功。如果匹配不成功： a：當(dāng)kK時(shí)，增加新的高斯模型； b：當(dāng)k=K時(shí)，用新高斯模型代替優(yōu)先級(jí)最小的模型。新的高斯模型，用當(dāng)前像素值作為新模型的均值，即，協(xié)方差為，權(quán)重為，其中為學(xué)習(xí)速率。（4） 未匹配模式的均值和方差

7、不變，對(duì)匹配模式的第i個(gè)高斯模型參數(shù)進(jìn)行更新： （5） 高斯模型參數(shù)更新完畢后，對(duì)每個(gè)像素點(diǎn)的K歌高斯模型按優(yōu)先級(jí)降序排序。取前B個(gè)高斯模型作為背景像素的最佳描述： （6） 繼續(xù)對(duì)與上述B個(gè)高斯模型進(jìn)行匹配檢驗(yàn)，如果與前B個(gè)高斯模型的任意一個(gè)匹配，則該像素點(diǎn)為背景點(diǎn)；否則為前景點(diǎn)。（7） 重復(fù)步驟（2）-（6），直到視頻結(jié)束。5、 GMM代碼實(shí)現(xiàn)#include#include#includeusing namespace cv;using namespace std;#define COMPONET 5 /混合高斯模型個(gè)數(shù)#define ALPHA 0.03 /學(xué)習(xí)率#define SD_I

8、NIT 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 *weight, *mean, *sigma, *u_diff, *rank;double p = ALPHA / (1 / (double)CO

9、MPONET);double rank_temp = 0;int rank_index_temp = 0;CvRNG state; /隨機(jī)生成狀態(tài)器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 = c

10、vCreateImage(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); /優(yōu)先級(jí)weight = (double*)malloc(siz

11、eof(double)*width*height*COMPONET); /權(quán)重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 fro

12、m mean /初始化均值、方差、權(quán)重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*width*COMPONET + j*COMPONET + k = SD_INIT;weighti*width*COMPONET + j*COMPONET + k = (double)1 / COMPONET;while (1)rank_in

13、dex = (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 imageDatai*width + j - meani*width*COMPONET + j*COMPONET + k);/update gaussian components fo

14、r 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_diffi*width*COMPONET + j*COMPONET + k imageDatai*width + j;sigmai*width*COMPONET + j*COMPONET + k = sqrt(1 - p)*(sigmai*width*COMPONET + j*COMPONET + k * sigmai*

15、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-imageDatai*width + j;sigmai*width*COMPONET + j*COMPONET + k

16、= 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 = (1 - ALPHA)*weighti*width*COMPONET + j*COMPONET + k; / wei

17、ght slighly decreases /sum_weight += weighti*width*COMPONET + j*COMPONET + k;/權(quán)重歸一化for (int k = 0; k COMPONET; k+)weighti*width*COMPONET + j*COMPONET + k = weighti*width*COMPONET + j*COMPONET + k / sum_weight;/獲取權(quán)重最小下標(biāo)double temp = weighti*width*COMPONET + j*COMPONET;int min_index = 0;backgi*width +

18、 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 imageDatai*width + j = (uchar)backgi*width + j;/if no components match, create new component if (mat

19、ch = 0)meani*width*COMPONET + j*COMPONET + min_index = (uchar)grayFrame-imageDatai*width + j;sigmai*width*COMPONET + j*COMPONET + min_index = SD_INIT;weighti*width*COMPONET + j*COMPONET + min_index = 1 / COMPONET;/計(jì)算優(yōu)先級(jí)for (int k = 0; k COMPONET; k+)rankk = weighti*width*COMPONET + j*COMPONET + k /

20、sigmai*width*COMPONET + j*COMPONET + k;/sort rank values for (int k = 1; k COMPONET; k+)for (int m = 0; m 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)