高斯混合模型视频背景建模的EM算法与Matlab 实现

1.问题描述

影像的背景前景分离.
输⼊为影像监控的1000 帧 (如下⽅图中左边所⽰), 要求输出是背景和前景 (如下⽅图中右边所⽰).

2.背景知识

观察待处理的监控影像，可以发现，前景主要是来来往往的行人，背景始终是摄像头对准的固定区域，可以认为是始终不变的，所以进行影像背景前景分离的一个思路是先对背景进行建模，获得对背景的模型描述后，前景也就不难获得了。
首先，我们做出假设：在整个视频中，每一帧的背景都是相同的。
严格的讲，视频每帧之间的背景不可能完全相同的，或多或少都会有一些像素点灰度发生变化，但是在我们的问题背景下，这种程度的差异是可以忽略的，认为每一帧的背景相同得到的效果是可以接受的。
基于这样的假设，我们可以对背景建模，并用每一帧减去同一个背景模型，来获得前景的信息。

2.1混合高斯模型

背景建模的方法,这里尝试经典的混合高斯模型。
混合高斯模型的思想源于单高斯模型。
在单高斯模型方法中，对于每一个像素，用一个高斯分布来描述其在不同时刻的灰度分布情况。
之所以可以用高斯分布描述，是基于假设：背景在图像序列中总是最经常被观测到。
然而，单一的高斯分布描述能力是很有限的，面对背景内容中的一些信息变化是不够敏感的，混合高斯模型应运而生。
一个高斯分布不够，就用多个高斯分布加权求和来描述一个像素点在不同时刻的灰度分布情况。

定义：高斯混合模型：具有如下形式的概率分布模型

P(y|θ)=∑k=1Kαkϕ(y|θ)

其中

αk是系数，∑Kk=1αk=1;

ϕ(y|θ)是高斯分布密度，θk=(μk,σ2k)

ϕ(y|θ)=12π√σkexp(−(y−μk)22σ2k)

称为第 k 个分模型。

使用上述公式描述的混合高斯模型理论上可以更准确地对视频背景进行建模。
确认使用的数学模型后，主要任务就是确认模型参数。
混合高斯模型中，对于任意一个分模型，有三个待确定参数：权值αj ,均值μ，方差σ2,一个混合高斯包含K个分模型，意味着需要确定3∗K个参数。
对于视频中尺寸为m∗n的帧图片而言，每一个像素点对应一个混合高斯模型的话，就需要确定m∗n∗K∗3个参数。
参数的选择对于模型的描述能力起着决定性作用，这里，我们选择使用EM算法迭代进行混合高斯模型的参数估计。

2.2EM算法

EM算法，期望极大算法，expectation maximization algorithm
EM 算法的每次迭代由两步组成：E步，求期望（expectation）；M 步，求极大（maximization）。
关于EM算法的具体介绍，参考李航老师的《统计学习方法》第9章。

高斯混合模型参数估计的EM算法
输入：观测数据y1,y2,...,yN,高斯混合模型
输出：高斯混合模型参数
(1)取参数的初始值开始迭代
(2)E步：依据当前模型参数，计算分模型k对观测数据 yj 的响应度

γ^jk=αkϕ(yj|θk)∑Kk=1αkϕ(yi|θ),j=1,2,...,N;k=1,2,...,K

(3)M步：计算新一轮迭代的模型参数

μ^k=∑Nj=1γ^jkyi∑Nj=1γ^jk,k=1,2,...,K

σ^2k=∑Nj=1γ^jk(yi−μk)2∑Nj=1γ^jk,k=1,2,...,K

α^k=∑Nj=1γ^jkN,k=1,2,...,K

(4)重复第(2)步和第(3)步，直到收敛

收敛的判断通过关于参数θ的对数极大似然函数Q函数进行。

Q(θ,θ(i))=∑k=1K{nklogαk+∑k=1Nγ^jk[log(12π−−√)−logσk−12σ2k(yj−μk)2]}

Q 函数描述的是对应高斯混合模型与数据点之间的接近程度，在EM算法的迭代过程中，计算每一次迭代得到的模型对应的Q函数值，当对数似然函数的变换小于预设阈值，也就是不再明显变化时，迭代结束。

3.Matlab实现

(1)Q函数的计算

function Q=qfunction(alpha,sigma,mu,D,K,fmax,i,gama)
%gama=gama(D,mu,alpha,sigma,k,i,fmax);
a = zeros(1,K);
b = zeros(1,K);
Q = zeros(1,fmax);
for j = 1:fmax                     %只取前fmax帧建模背景
    for m = 1:K
        a(m)=(sum(gama(:,m)))*log(alpha(m));
        b(m)=gama(j,m)*(log(1/sqrt(2*pi))-log(sigma(m))-((D(i,j)-mu(m))^2)/(2*(sigma(m))^2));   
    end
    Q(j) = sum(a+b);
end

(2)响应度γ的计算

function gama=fgama(D,mu,alpha,sigma,k,i,fmax)
sgama =zeros(1,k);
ogama = zeros(1,k);
gama = zeros(fmax,k);
for e = 1:fmax
    for c = 1:k
        gauss = 1/sqrt(2*pi)*exp((-(D(i,e)-mu(c))^2)/(2*sigma(c))^2);
        sgama(c)=alpha(c)*gauss;   
    end
    bgama = sum(sgama);
    for d = 1:k
    ogama(d)=sgama(d)/bgama;
    end
    gama(e,:)=ogama;
end

(3)视频帧的读入

path(path,genpath(pwd)); 
files = dir('hall-partial*.bmp');
 if numel(files)<10; 
     error('not enough frames, at least 10 frames');
 end;
 H = imread(files(1).name);
 [row,colum,chan] = size(H);
 imsize = [row,colum];   D = [];
 if chan>1 % if multi channel, only take one channel
   RGBchan = 0; % we can use second (green) channel, but we can also use other channel
    for k = 1:min(200,numel(files))
      H = imread(files(k).name);
      if RGBchan
      H = H(:,:,RGBchan);
      else
      H = rgb2gray(H);    
      end
      D(:,k)=H(:);
    end
 else
   for k = 1:min(200,numel(files))
     H = imread(files(k).name);
     H = H(:,:,1);
     D(:,k)=H(:);
   end   
 end

(4)背景建模主程序

%-----------------混合高斯模型背景建模--------------------%
 %--------------------------------------------------------%
 pn = row * colum;           %像素数量 
 K = 3;                     %分模型个数为K个
 fmax = 20;
 mu = zeros(K,1);            %均值
 sigma = zeros(K,1);         %方差
 alpha = zeros(K,1);         %权值
 BM = zeros(pn,K);        %均值矩阵
 BM2 = zeros(pn,1);
 t = 1;                      %阈值

 %gama = zeros(fmax,k);
 for i = 1:pn
     for k = 1:K
         mu(k) = D(i,1);
         sigma(k) = 20+(k-1)*1 ;
     end
     for k = 1:K
         alpha(k) = 1-(sigma(k)/sum(sigma));         %初始化起始参数
     end
     gama = fgama(D,mu,alpha,sigma,k,i,fmax);
     newQ = qfunction(alpha,sigma,mu,D,K,fmax,i,gama);
     oldQ = zeros(size(newQ));
     while (abs(sum(oldQ-newQ)) >= t*fmax)
         for k = 1:K
             sg = sum(gama(:,k));
             mu(k) = (sum(gama(:,k).*D(i,1:fmax)'))/sg;
             sigma(k) = (sum(gama(:,k).*((D(i,1:fmax)'-repmat(mu(k),fmax,1)).^2)))/sg;
             alpha(k)=sg/fmax;
         end
         gama = fgama(D,mu,alpha,sigma,k,i,fmax);
         oldQ = newQ;
         newQ = qfunction(alpha,sigma,mu,D,K,fmax,i,gama);
     end
     muu = sum(alpha.*mu);
     sigmau = sum(alpha.*mu);
     BM(i,:)=mu';
     BM2(i,:)=muu;    
 end

 BP = repmat(BM2(:,1),1,200);
 FP = D - BP;
 FP(abs(FP)<1.5*sigmau)=0;
 mat2im(BP,imsize);
 if ~exist('frontpeople.avi','file')
     mat2movG(FP,imsize,1);
 end

(5)结果输出部分的代码与主题关系不大，就不贴出了。

4.运行结果

分离得到的背景图片如下

分离得到的前景视频缩略图如下

5.分析总结

首先，程序得到的背景前景分离结果是不错的，效果满足需求，说明混合高斯模型的理论效果得到了证明，

不足的地方主要在于程序的运行效率。由于混合高斯模型待确定参数非常多，有数万个，而每个参数的确定都需要经过EM算法的数次迭代，每次迭代过程中的计算量也非常大，几个因素导致程序的运行时间非常长。在调整过初值和终止阈值等参数后，也需要一分钟左右的运行时间，这样的运行效率是无法接受的。

显然，EM算法不是求解混合高斯模型的最高效算法，但EM算法本身具有广泛的应用范围，所以这里使用EM算法更多地是为了熟悉学习算法流程。
混合高斯模型的参数估计还有经典的实时更新的算法。计算效率会高很多，由于没有自己动手实现过，就不赘述了。

撰写用时：2小时
参考文献 ：《统计学习方法》 李航 第9章
Written with StackEdit.