混合高斯模型用于动态背景建模

ref: Improved adaptive Gausian mixture model for background subtraction,Zoran Zivkovic,2004

摘要

背景提取是一个常见的计算机视觉问题，我们分析了一般像素级方法，提出了一种有效的自适应高斯混合概率密度方法，使用回归方程更新模型参数，同时为每个像素选取合适的模型数量。

简介

一种自上而下的方法就是为每个像素点建立一个概率密度函数。对于新图像中一个像素，如果它的值能够被概率密度函数很好地描述，那么它属于背景。对于一个固定的场景，最简单的背景模型可以是尚不存在侵入物体时的场景图像。下面要做的就是为每个像素估计合适的方差，不同像素的方差一般不同。文献[1] 中就使用了单个高斯分布模型。然而像素值的分布一般比较复杂，所以需要更加复杂的描述模型。文献[2] 中提出的使用混合高斯分布模型GMM用于背景提取。文献[3]提出了现在常用的更新 GMM 模型参数的方法，文献[10] 提出了更为复杂的更新方法。这些GMM算法使用固定的模型个数。我们对此基于文献[12] 提出了一种改进的GMM算法，每个像素的高斯分布参数和高斯模型个数都会自适应更新。通过在线的更新每个像素的高斯模型个数，算法能够自适应地满足场景需求。

混合高斯分布模型

实际场景中光照会逐渐变化（室外的日照或者天气原因）或者突然变化（室内灯光的开闭），场景中也存在物体的突然出现或消失。为了适应这些变化，我们可以通过增加新样本替换旧样本来更新训练集。我们选择一个合适的时间窗口 T，所以时间 t 时的样本集 (chi_T = {x^t,...,x^{(t-T)}}) 。对于每个新样本我们更新 (chi _T) 并重新估计 (hat{p}(vec{x}|chi,BG+FG))。我们使用具有 M 个分布的混合高斯模型，如下式所示：

[hat{p}(vec{x}|chi_T,BG+FG)=sum^M_{m=1} hat{pi}_m aleph(vec{x};hat{vec{mu}}_m,hat{sigma}^2_m I) (3) ]

其中 (hat{mu}_1,...,hat{mu}_M) 是高斯分布均值，(hat{sigma}_1,...,hat{sigma}_M) 是高斯分布方差估计值。协方差矩阵为对角矩阵（其中假设 RGB 三通道数据是互相独立的且具有相同的方差。虽然这不符合实际情况，但我们通过牺牲一定精度来避免复杂的矩阵运算）。(hat{pi}_m) 表示的权重非负且和为1。时间 t 时的像素值为 (vec{x}^t)，模型更新方程如下所示：

[egin{cases} hat{pi}_m leftarrow hat{pi}_m + alpha(o^t _m - hat{pi}_m) \hat{vec{mu}}_m leftarrow hat{vec{mu}}_m + o_m ^t (alpha / hat{pi}_m )vec{delta}_m \hat{sigma}^2_m leftarrow hat{sigma}_m ^2 + o_m^t (alpha /hat{pi}_m)(vec{delta}^T _m vec{delta}_m - hat{sigma}^2_m) & (4)-(6) end{cases} ]

其中 (vec{delta}=vec{x}^t - hat{vec{mu}}_m)。这里我们没有使用之前提到的 T，而是使用 (alpha) 表示指数级下降的包络曲线，用以限制旧数据的影响，我们可以近似理解为 (alpha = 1/T)。

由于图像中每一个像素都建立了混合高斯模型，使用 EM 算法进行精确求解时会很复杂，这里使用了一种类似于 K 均值的方法来近似求解。

对于一个新样本，它最符合的分布模型的归属因子 (o^t _m) 为 1，其他分布模型的归属因子为 0。匹配的条件定义为新像素值位于高斯分布的 3 倍标准差范围内。如果像素值不符合所有高斯分布，会建立一个新的高斯分布，其中 (hat{pi}_{M+1}=alpha,hat{vec{mu}}_m = vec{x}^t, vec{sigma}_{M+1} = sigma_0)，(sigma_0) 是一个合适的初始值。如果分布中模型个数达到最大值，我们会把 (vec{pi}_m) 最小的分布舍弃，同时引入以当前像素值为均值，权重为一个较低的初始值，方差为一个较高的初始值的高斯分布。

这种方法的一点好处在于当某个区域被计算归为背景时，它不会破坏当前区域的高斯分布。之前的背景高斯分布会继续存在于混合模型中，直至它成为按权重降序排列的第 k 个高斯模型同时一个不与当前混合模型匹配的新的像素值被观测到。所以如果一个物体静止了足够长的时间使得它被计算归为背景区域后开始运动，那么之前癿背景分布模型会仍然存在，它的 (mu) 和 (sigma) 保持不变，权值会降低，同时会很快再次被归类为背景。

背景模型估计

每个像素的混合高斯分布参数会一直在变，我们想要确定其中哪一个高斯分布更加符合当前的背景。经验主义上来说，我们倾向于选择权重更大的同时方差更小的高斯分布。

为什么有这样的结论，我们做如下解释。当一个物体固定存在时，会产生一个高斯分布的权重累积增加和方差累积减小。作为对比，当一个新的物体覆盖了原有区域时，一般它不会匹配当前的混合高斯分布模型，所以要么产生一个权重很小方差很大的新的高斯分布，或者会使当前的一个高斯分布的方差变大。同时一般来说一个运动物体产生的分布的方差要比静止物体的要大，直到物体停止运动。为了精确描述，我们需要定量的方法来确定当前混合模型中哪部分分布能够最好的表达当前背景。

该算法中使用一种在线聚类方法，一般来说新出现的前景物体会由新增加的几个权重比较小的聚类（高斯分布）表示，所以我们可以近似使用前B 个聚类表达背景模型：

[hat{p}(vec{x}|chi_T,BG) sim sum^B_{m=1} hat{pi}_m aleph(vec{x};hat{vec{mu}}_m,hat{sigma}^2_m I) (7) ]

如果通过 (hat{pi}_m) 降序排列，可以得到

[B={argmin_b} (sum_{m=1} ^b hat{pi}_m > 1-c_f) (8) ]

这里 (c_f) 是像素值被归为前景物体且不影响背景模型的最大比率。举例说明如下，如果一个新的物体出现在场景中且持续了一定时间，它会产生一个新的比较稳定的分布。由于物体不满足以前的分布，新分布的权值 (pi _{B+1}) 会逐渐增加。如果物体出现的时间足够长，分布权值会增加到比 (c_f) 还要大，它就被看做是背景模型的一部分。我们观察公式（4）后可以知道物体需要出现 (log(1-c_f)/log(1-alpha)) 帧后才会满足上述情况。

模型中分布个数的选择

权值 (pi_m) 表述的是数据有多大可能属于 GMM 的第 m 个高斯分布。我们可以认为它定义了数据来自第 m 个分布的概率，所以 (pi _m) 序列定义了一个隐式的多项式分布。现在我们假设有 t 个数据采样值并且都属于 GMM 中的分布。我们同时假设属于第 m 个分布的样本数为 (n_m = sum^t _{i=1}o^i_m)，(o^i _m) 定义参考上文。对于 (n_m) 假设的多项式分布可以得到似然函数 (L=prod^M_{m=1}pi^{n_m}_m) （这里涉及到了最大似然估计相关的内容，基本思想就是分布模型的参数应该使得样本出现的概率最大）。同时存在约束条件就是权值的和为 1，我们通过引入拉格朗日惩罚因子 (lambda) 来使用该约束条件。最大似然估计计算得到 ({partial over partial hat{pi}_m}(logL+lambda(sum^M_{m=1}hat{pi}_m-1))=0) （就是导数为 0 的地方似然函数取值最大，引入 log 的原因是方便求导）。通过对 (lambda) 处理后得到

[hat{pi}^t_m = frac{n_m}{t} = frac{1}{t}sum^t_{i=1}o^i_m (9) ]

这部分推到过程如下：

[{partial over partial hat{pi}_m}(logL+lambda(sum^M_{m=1}hat{pi}_m-1))=0 ]

得到

[frac{n_m}{hat{pi}_m}+lambda = 0 ightarrow hat{pi}_m =frac{-n_m}{lambda}(1le m le M) ]

[frac{hat{pi}_1 + ... + hat{pi}_m}{n_1 + ... + n_m} = frac{1}{t} = frac{1}{-lambda} ightarrow lambda = -t ]

(hat{pi}^t_m) 可以写成递归的形式，得到

[hat{pi}_m^t = hat{pi}_m^{t-1} + 1/t(o^t _m - hat{pi}_m^{t-1}) (10) ]

如果我们将 1/t 设定为一个定值（固定了新采样数据的影响），譬如就取为 (alpha = 1/T) ，那我们就得到了更新方程（4）。这种新采样数据固定的影响表明我们更多以新采样值为参考更新权值，以前的采样值影响权值会像上文提到的那样指数级衰减。

这里也补充方程（5） 的推导过程。

（后面的内容涉及到狄利克雷分布，在最大似然估计过程中加入对参数分布先验概率密度函数的考虑，相关内容可参考狄利克雷分布讲解 ，如何理解 Beta 分布和 Dirichlet 分布 ，先验概率、后验概率以及共轭先验 ）

基于多项式分布的研究结果，我们引入它的共轭先验，狄利克雷分布 (P=prod^M_{m=1}pi ^{c_m}_m) 。因子 (c_m) 对于多项式分布来说，从最大后验估计意义上讲就是第 m 个分布的先验，即属于第 m 个分布的先验样本数量。就像在文献[12] 中的那样，我们取 (c_m) 为负数，即 (c_m = -c) 。这表明我们只有在存在足够多样本支撑某个分布的存在时才会接受它（个人理解这个地方取负值更多的是为了求导时的符号变化）。这种先验也与文献[12] 中用于选择合适模型的采样值的最小信息长度相关。最大后验估计方法中也包含求导过程，即

[frac{partial}{partial hat{pi}_m}(logL+logP+lambda(sum^M_{m=1}hat{pi}_m-1)) = 0 ]

其中 (P=prod^M_{m=1}pi ^{-c}_m) ，计算后我们得到：

[hat{pi}^t_m=frac{1}{K}(sum^t_{i=1}o^i_m -c) (11) ]

其中 (K=sum_m(sum_to^t_m -c) = t-Mc) 。我们可以重写公式（11）为

[hat{pi}^t_m = frac{hat{Pi}_m-c/t}{1-Mc/t} (12) ]

其中 (hat{Pi}_m=frac{1}{t}sum^t_{i=1}o^i_m) ，即公式（9）根据最大似然估计得到的结果。公式（12）中通过 (c/t) 引入了先验估计的影响（这里应该是 (c_m=-c) 的原因），当样本集合更大时该影响变小。不过如果一个小的偏置可接受的话，我们可以把 (c/t) 固定为一个常数，即 (c/t=c_T=c/T) ，这里 T 取较大的值。这表明对于有 T 个样本的集合，偏置影响将是一个固定量。然后把这个带入公式（11）写为递归的形式如下：

[hat{pi}^t_m = hat{pi}^{t-1}_m + 1/t(frac{o^t_m}{1-Mc_T}-hat{pi}^{t-1}_m)-1/t frac{c_T}{1-Mc_T} (13) ]

一般情况下 M 取值不大且 (c_T) 也很小，我们可以得到 (1-Mc_T approx 1) 。然后把 (1/t) 替换成 (alpha) ，我们得到最终的更新方程如下：

[hat{pi}_m leftarrow hat{pi}_m + alpha(o^t _m - hat{pi}_m)-alpha c_T (14) ]

与公式（4）相比，通过最大后验估计引入了一个减常量，这是本文的最大创新点，之后根据权值为负时会把这个分布舍弃掉，实现分布数量的自适应性）

该公式将替代公式（4）使用，每次更新后我们需要将权值 (pi_m) 归一化。我们将 GMM 初始化为一个高斯分布，该分布的均值是像素值，之后参考文献[3] 增加新的高斯分布。负数权值的狄利克雷分布将使与该像素不匹配的分布权值减小，当其权值为负时我们舍弃该分布，这也保证了分布权值式中非负。对于 (alpha = 1/T) ，我们设定至少要有 (c=0.01 imes T) 个采样值符合一个分布，所以这里 (c_t=0.01) 。

算法实现

当然 opencv 已经打包了 MOG，感兴趣的可以学习源码。

//------------------MOG_BGS.h------------------
#pragma once  
#include <iostream>  
#include "opencv2/opencv.hpp"  

using namespace cv;
using namespace std;

//定义gmm模型用到的变量  
#define GMM_MAX_COMPONT 6        //高斯模型个数  
#define GMM_LEARN_ALPHA 0.005    //学习速率
#define GMM_THRESHOD_SUMW 0.7    //前景背景分割阈值

class MOG_BGS
{
public:
	MOG_BGS(void);
	~MOG_BGS(void);

	void init(const Mat _image);
	void processFirstFrame(const Mat _image);
	void trainGMM(const Mat _image);
	void getFitNum(const Mat _image);
	void testGMM(const Mat _image);
	Mat getMask(void){ return m_mask; };

private:
	Mat m_weight[GMM_MAX_COMPONT];  //权值  
	Mat m_mean[GMM_MAX_COMPONT];    //均值  
	Mat m_sigma[GMM_MAX_COMPONT];   //方差  

	Mat m_mask;                     //前景掩模
	Mat m_fit_num;
};


//------------------MOG_BGS.cpp------------------
#include "MOG_BGS.h"

MOG_BGS::MOG_BGS(void)
{

}

MOG_BGS::~MOG_BGS(void)
{

}

// 全部初始化为0
void MOG_BGS::init(const Mat _image)
{
	/****initialization the three parameters ****/
	for (int i = 0; i < GMM_MAX_COMPONT; i++)
	{
		m_weight[i] = Mat::zeros(_image.size(), CV_32FC1);
		m_mean[i] = Mat::zeros(_image.size(), CV_8UC1);
		m_sigma[i] = Mat::zeros(_image.size(), CV_32FC1);
	}
	m_mask = Mat::zeros(_image.size(), CV_8UC1);
	m_fit_num = Mat::ones(_image.size(), CV_8UC1);
}

//gmm第一帧初始化函数实现
//捕获到第一帧时对高斯分布进行初始化．主要包括对每个高斯分布的权值、期望和方差赋初值．
//其中第一个高斯分布的权值为1，期望为第一个像素数据．其余高斯分布权值为0，期望为0．
//每个高斯分布都被赋予适当的相等的初始方差 15
void MOG_BGS::processFirstFrame(const Mat _image)
{
	for (int i = 0; i < GMM_MAX_COMPONT; i++)
	{
		if (i == 0)
		{
			m_weight[i].setTo(1.0);
			_image.copyTo(m_mean[i]);
			m_sigma[i].setTo(15.0);
		}
		else
		{
			m_weight[i].setTo(0.0);
			m_mean[i].setTo(0);
			m_sigma[i].setTo(15.0);
		}
	}
}

// 通过新的帧来训练GMM
void MOG_BGS::trainGMM(const Mat _image)
{
	for (int i = 0; i < _image.rows; i++)
	{
		for (int j = 0; j < _image.cols; j++)
		{
			int num_fit = 0;

			/**************************** Update parameters Start ******************************************/
			for (int k = 0; k < GMM_MAX_COMPONT; k++)
			{
				int delm = abs(_image.at<uchar>(i, j) - m_mean[k].at<uchar>(i, j));
				long dist = delm * delm;
				// 判断是否匹配：采样值与高斯分布的均值的距离小于3倍标准差（9倍方差）
				if (delm < 9.0 * m_sigma[k].at<float>(i, j))
				{
					// 如果匹配
					/****update the weight****/
					m_weight[k].at<float>(i, j) += GMM_LEARN_ALPHA * (1 - m_weight[k].at<float>(i, j));

					/****update the average****/
					m_mean[k].at<uchar>(i, j) += (GMM_LEARN_ALPHA / m_weight[k].at<uchar>(i, j)) * delm;

					/****update the variance****/
					m_sigma[k].at<float>(i, j) += (GMM_LEARN_ALPHA / m_weight[k].at<float>(i, j)) * (dist - m_sigma[k].at<float>(i, j));
				}
				else
				{
					// 如果不匹配。则该该高斯模型的权值变小
					m_weight[k].at<float>(i, j) += GMM_LEARN_ALPHA * (0 - m_weight[k].at<float>(i, j));
					num_fit++; // 不匹配的模型个数
				}
			}
			/**************************** Update parameters End ******************************************/


			/*********************** Sort Gaussian component by 'weight / sigma' Start ****************************/
			//对gmm各个高斯进行排序,从大到小排序,排序依据为 weight / sigma
			for (int kk = 0; kk < GMM_MAX_COMPONT; kk++)
			{
				for (int rr = kk; rr< GMM_MAX_COMPONT; rr++)
				{
					if (m_weight[rr].at<float>(i, j) / m_sigma[rr].at<float>(i, j) > m_weight[kk].at<float>(i, j) / m_sigma[kk].at<float>(i, j))
					{
						//权值交换
						float temp_weight = m_weight[rr].at<float>(i, j);
						m_weight[rr].at<float>(i, j) = m_weight[kk].at<float>(i, j);
						m_weight[kk].at<float>(i, j) = temp_weight;

						//均值交换
						uchar temp_mean = m_mean[rr].at<uchar>(i, j);
						m_mean[rr].at<uchar>(i, j) = m_mean[kk].at<uchar>(i, j);
						m_mean[kk].at<uchar>(i, j) = temp_mean;

						//方差交换
						float temp_sigma = m_sigma[rr].at<float>(i, j);
						m_sigma[rr].at<float>(i, j) = m_sigma[kk].at<float>(i, j);
						m_sigma[kk].at<float>(i, j) = temp_sigma;
					}
				}
			}
			/*********************** Sort Gaussian model by 'weight / sigma' End ****************************/


			/*********************** Create new Gaussian component Start ****************************/
			if (num_fit == GMM_MAX_COMPONT && 1.0e-6 > m_weight[GMM_MAX_COMPONT - 1].at<float>(i, j))
			{
				//if there is no exit component fit,then start a new component
				//当有新值出现的时候，若目前分布个数小于M，新添一个分布，以新采样值作为均值，并赋予较大方差和较小权值
				for (int k = 0; k < GMM_MAX_COMPONT; k++)
				{
					if (1.0e-6 > m_weight[k].at<float>(i, j))
					{
						m_weight[k].at<float>(i, j) = GMM_LEARN_ALPHA;
						m_mean[k].at<uchar>(i, j) = _image.at<uchar>(i, j);
						m_sigma[k].at<float>(i, j) = 15.0;

						break; //找到第一个权值不为0的即可
					}
				}
			}
			else if (num_fit == GMM_MAX_COMPONT && m_weight[GMM_MAX_COMPONT - 1].at<float>(i, j) != 0)
			{
				//如果GMM_MAX_COMPONT都曾经赋值过，则用新来的高斯代替权值最弱的高斯,权值不变
				m_mean[GMM_MAX_COMPONT - 1].at<uchar>(i, j) = _image.at<uchar>(i, j);
				m_sigma[GMM_MAX_COMPONT - 1].at<float>(i, j) = 15.0;

				for (int q = 0; q < GMM_MAX_COMPONT; q++)
				{
					//对权值进行更新，保持权值和为1  
					m_weight[q].at<float>(i, j) /= (1 - GMM_LEARN_ALPHA);
				}
			}
			/*********************** Create new Gaussian component End ****************************/
		}
	}
}

//对输入图像每个像素gmm选择合适的高斯分量个数
//排序后最有可能是背景分布的排在最前面，较小可能的短暂的分布趋向于末端．我们将排序后的前fit_num个分布选为背景模型;
//在排过序的分布中，累积概率超过GMM_THRESHOD_SUMW的前fit_num个分布被当作背景模型，剩余的其它分布被当作前景模型．
void MOG_BGS::getFitNum(const Mat _image)
{
	for (int i = 0; i < _image.rows; i++)
	{
		for (int j = 0; j < _image.cols; j++)
		{
			float sum_w = 0.0;	//重新赋值为0，给下一个像素做累积
			for (uchar k = 0; k < GMM_MAX_COMPONT; k++)
			{
				sum_w += m_weight[k].at<float>(i, j);
				if (sum_w >= GMM_THRESHOD_SUMW)
				{
					m_fit_num.at<uchar>(i, j) = k + 1;
					break;
				}
			}
		}
	}
}

//gmm测试函数的实现
void MOG_BGS::testGMM(const Mat _image)
{
	for (int i = 0; i < _image.rows; i++)
	{
		for (int j = 0; j < _image.cols; j++)
		{
			int k = 0;
			for (; k < m_fit_num.at<uchar>(i, j); k++)
			{
				if (abs(_image.at<uchar>(i, j) - m_mean[k].at<uchar>(i, j)) < 3 * sqrt(m_sigma[k].at<float>(i, j)))
				{
					m_mask.at<uchar>(i, j) = 0;
					break;
				}
			}
			if (k == m_fit_num.at<uchar>(i, j))
				m_mask.at<uchar>(i, j) = 255;
		}
	}
}