混和高斯模型的推導和實現(xiàn)..docx

1、基于 GMM的運動目標檢測方法研究一、GMM數(shù)學公式推導1、預備知識：（1）設(shè)離散型隨機變量X 的分布率為：P Xakpk ,k1,2則稱 E Xak p k 為 X 的數(shù)學期望或均值k（ 2）設(shè)連續(xù)型隨機變量X 的概率密度函數(shù)（ PDF）為 f(x)其數(shù)學期望定義為： E Xxf x dx（ 3） D XE XE X2稱為隨機變量 x 的方差，D X稱為 X的標準差（4）正態(tài)分布： X N,2概率密度函數(shù)為： p x1e2x222（ 5）設(shè) (x,y) 為二維隨機變量，EXE XYE Y若存在，則稱其為 X 和 Y 的協(xié)方差，記為 cov(x,y)cov X ,YE XE XYE YE XY

2、2、單高斯模型： SGM（也就是多維正態(tài)分布）其概率密度函數(shù)PDF定義如下：11xT C 1 xN x; ,Ce 22n C其中，x 是維數(shù)為 n 的樣本向量（列向量）， 是期望， C是協(xié)方差矩陣， |C|表示 C 的行列式， C 1 表示 C 的逆矩陣，xT 表示 x的轉(zhuǎn)置。3、混合高斯模型： GMM設(shè)想有 m 個類：1，2，3， ，m，每類均服從正態(tài)分布。各分布的中心點（均值）分別為：1， 2， 3， ， m方差分別為：1， 2， 3， ， m每一類在所有的類中所占的比例為P1 , P2 , P3 , Pmm其中Pi1。i1同時，已知個觀察點：。其中，用大寫P 表示概率，用小寫p 表示概率

3、密度。則依此構(gòu)想，可得概率密度函數(shù)為：p x N1,1P1N2,2P2N m , m P mmP1T1ixCxe 2i 12d C其中 d 是維數(shù)， | | 是行列式但是在利用 GMM進行目標檢測時， 這些模型的參數(shù)可能已知， 也可能不知道， 當參數(shù)已知時， 可以直接利用 GMM進行目標檢測， 在未知的情況下， 需要對參數(shù)進行估計。 對參數(shù)估計時，還要考慮樣本分類是否已知。（ 1）樣本已知：最大似然估計：可以直接采用 MLE（最大似然估計）進行參數(shù)估計：未知量為集合：1， ， m，C1， , Cm，P1 ， , P mn將衡量概率密度函數(shù)優(yōu)劣的標準寫出： p x |P xk |k 1即為：nm

4、p x |k 1 i1Pi1xkTC1xkie 2i2d | C |只要定出該標準的最大值位置，就可以求出最優(yōu)的待定參數(shù)。為了求出這個最nln px|ln kP xk |1nln Pxk |大值的位置，就需用導數(shù)求極點，具體求解過程于下：k1mnN xk, i P iln i 1k1mln px|nnln iN xk , iPi1k 1k 1m1m求導：mN xk ,iPi k 1N xk , i Pii 1i 1m1mPi1 xkiT C -1 xk ie2pxk|di12|C |k 1然后再分別對各個參數(shù)求導：求參數(shù)i：對感興趣，求偏導數(shù)有：對感興趣，接下來的求導比較復雜，在此就沒有繼續(xù)推

5、導。（ 2）樣本未知：EM 估計，算法流程：初始化：方案 1：協(xié)方差矩陣 C j 0 設(shè)為單位矩陣，每個模型比例的先驗概率設(shè)為j 01，M均值j 0 為隨機數(shù)。方案 2：有 K 均值（ K-means)聚類算法對樣本進行聚類，利用各類的均值作為j 0 ，并計算 C j 0 ，j 0 去各類樣本占總數(shù)的比例。估計步驟（ E-step ）：令 j 的后驗概率為：j N jxi |ijM,1 i n,1 j Mk Nkxi |k 1最大化步驟（ M-step ）：nij更新權(quán)值：i1jnnxiij更新均值：i1jniji1nTij xii xiii1更新方差矩陣： C jniji1收斂條件：不斷地迭

6、代步驟和， 重復更新上面的三個值， 直到 p X |p X | |，其中為更新參數(shù)后計算的值， 即前后兩次迭代得到的結(jié)果變化小于一定程度則終止迭代，通常10 -5二、 GMM發(fā)展歷史及現(xiàn)狀背景建模方法有很多種，如中值法、均值法、卡爾曼濾波器模型、碼本背景模型等，其中混合高斯模型是最經(jīng)典的算法。GMM最早是由 CHris Stauffer等在 1 中提出的，該方法是按照高斯分布對每個像素建立模型，并通過基于回歸濾波的在線 EM 近似方法對模型參數(shù)進行更新，它能魯棒地克服光照變化、樹枝搖動等造成的影響，但該方法也存在一些問題： 1）該方法對運動物體在場景中停止不動或者長時間停止時檢測失效， 而且?guī)?/p>

7、有初始學習速度慢， 在線更新費時、計算量大； 2）無法完整準確地檢測大并且運動緩慢的運動目標，運動目標的像素點不集中， 只能檢測到運動目標的部分輪廓， 無法提取出目標對象的完整區(qū)域； 3）無法將背景顯露區(qū)域與運動目標區(qū)域很好地區(qū)分開； 4）當運動目標由靜止緩慢轉(zhuǎn)化為運動時，易將背景顯露區(qū)檢測為前景，出現(xiàn)“影子”現(xiàn)象。三、GMM缺點及改進方法針對上述問題，一些科學研究者又在GMM算法的基礎(chǔ)上做了很多的改進： 張、白等人 2 引入分塊思想， 把圖像分為 L*L 塊；黃、胡等人 3 也引入了分塊的思想，但是他們的分塊理念是以當前像素點的8 鄰域作為一塊；華、劉4 把GMM與改進的幀差法( 相鄰兩幀圖

8、像對應(yīng)像素點8 鄰域像素值相減之和) 相結(jié)合，提高了計算效率； Suo 等人 5 是將混合高斯模型中的模型個數(shù)采改進為自適應(yīng)的；劉等人6融合幀間差分法，檢測背景顯露區(qū)域和運動區(qū)域，很好的解決了問題4。除此之外，還有基于紋理的混合高斯模型。四、GMM算法流程（ 1）用第一幀圖像對高斯混合模型進行初始化0,0x yI x y0x, ystd _ init02x, ystd _ init std _ initw01M一般模型的個數(shù) M為 3-6個，其中 std_init設(shè)置為 20（ 2）對于 t 時刻的像素 I t x, y ，分別與已經(jīng)存在的M個高斯模型依次進行匹配：I t x, yi ,t 1

9、 ( x, y) | 2.5 i ,t 1（ 3）如果滿足匹配條件，則該像素值與高斯模型匹配成功。如果匹配不成功：a ：當 kK 時，增加新的高斯模型；b：當k=K 時，用新高斯模型代替優(yōu)先級i最小的模型。新的高斯模型，用當前像素值作為新模型的均值，即iI x, y，協(xié)方差為istd _ init，權(quán)重為wi，其中 為學習速率。（ 4）未匹配模式的均值和方差不變，對匹配模式的第 i 個高斯模型參數(shù)進行更新：i ,t1i,t 1Itx y,212I tx, yi ,t 12i,ti,t 1wi,t1wi ,t 1（ 5）高斯模型參數(shù)更新完畢后， 對每個像素點的K 歌高斯模型按優(yōu)先級i降序排序。取

10、前 B 個高斯模型作為背景像素的最佳描述：MBkk 1iT（ 6）繼續(xù)對 I t x, y 與上述 B 個高斯模型進行匹配檢驗，如果I t x, y 與前 B 個高斯模型的任意一個匹配，則該像素點為背景點；否則為前景點。（ 7）重復步驟（ 2）- （6），直到視頻結(jié)束。五、GMM代碼實現(xiàn)#include#include#includeusing namespace cv;using namespace std;#define COMPONET 5/混合高斯模型個數(shù)#define ALPHA 0.03 /學習率#define SD_INIT 6/方差初始值#define THRESHOLD 0.

11、25/前景所占比例#define D 2.5int 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)COMPONET);double rank_temp = 0;int r

12、ank_index_temp = 0;CvRNG state; /隨機生成狀態(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 = cvCreateImage(CvSize(frame-width, f

13、rame-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)先級weight=(double*)malloc(sizeof(double)*width*height*COMPONET); /權(quán)重m

14、ean=(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/ 初始化均值、方差、權(quán)重for (int i = 0; i height; i+)for (i

15、nt j = 0; j width; j+)for (int k = 0; k COMPONET; k+)meani*width*COMPONET +cvRandReal(&state) * 255;sigmai*width*COMPONET+j*COMPONET + j*COMPONET +kk=SD_INIT;weighti*width*COMPONET+j*COMPONET+k=(double)1 / COMPONET;while (1)rank_index = (int *)malloc(sizeof(int)*COMPONET); cvCvtColor(frame, grayFram

16、e, 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 for each pixel for (int i = 0; i height; i+)for (int j = 0; j width; j

17、+)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*width*COMPONET+j*COMPONET+k)+p*(pow(uchar)grayFrame-imageDatai*width +

18、 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 =sqrt(1 - ALPHA)*(sigmai*width*COMPONET + j*COMPONET + k *sigmai*width*COMPONET+j*C

19、OMPONET+k)ALPHA*(pow(uchar)grayFrame-imageDatai*width+j+-meani*width*COMPONET + j*COMPONET + k, 2);/else/ weighti*width*COMPONET (1 - ALPHA)*weighti*width*COMPONET + j*COMPONET+ k j*COMPONET+ k;=/ weight slighly decreases/sum_weight += weighti*width*COMPONET + j*COMPONET + k;/ 權(quán)重歸一化for (int k = 0; k

20、 COMPONET; k+)weighti*width*COMPONET + j*COMPONET+ k = weighti*width*COMPONET + j*COMPONET + k / sum_weight;/ 獲取權(quán)重最小下標double temp = weighti*width*COMPONET + j*COMPONET;int min_index = 0;backgi*width + j = 0;for (int k = 0; k COMPONET; k+)backgi*width + j = backgi*width meani*width*COMPONET + j*COMPO

21、NET +kj+*weighti*width*COMPONET + j*COMPONET + k;if (weighti*width*COMPONET + j*COMPONET + kimageDatai*width+j=(uchar)backgi*width + j;/if no components match, create new componentif (match = 0)meani*width*COMPONET+ j*COMPONET+ min_index= (uchar)grayFrame-imageDatai*width + j; sigmai*width*COMPONET

22、+ j*COMPONET +min_index = SD_INIT;weighti*width*COMPONET + j*COMPONET + min_index = 1 / COMPONET;/ 計算優(yōu)先級for (int k = 0; k COMPONET; k+)rankk= weighti*width*COMPONET + j*COMPONET+ k / sigmai*width*COMPONET + j*COMPONET + k;/sort rank valuesfor (int k = 1; k COMPONET; k+)for (int m = 0; m rankm)/swap max valuesrank_temp = rankm;rankm = rankk;rankk = rank_temp;/swap max index valuesrank_index_temp = rank_indexm;rank_indexm = rank_indexk;rank_index