【機器視覺】7. 目標跟蹤：背景減除、粒子濾波、meanshift

05-21

來自專欄 Let Machine See

目標跟蹤：背景減除、粒子濾波、meanshift

我們已經介紹了目標跟蹤中最基本的光流演算法

Calvin Shi：【機器視覺】3. 目標跟蹤：光流法?

zhuanlan.zhihu.com

以及卡爾曼濾波

Calvin Shi：【機器視覺】4. 卡爾曼濾波?

zhuanlan.zhihu.com

這篇文章接著介紹一些目標跟蹤中常用的演算法。

基於中值濾波的背景減除

背景減除（background subscription）目的在於區分場景背景和場景運動，從而進行運動檢測。

中值濾波器是順序統計濾波器的一種，即用該像素的相鄰像素的灰度中值來代替該像素的值，傳統中值濾波(Standard Median Filter)通過快速排序尋找像素點集合的中值，該方法對濾波窗口形狀和像素灰度級沒有限制，具有更加廣泛的運用範圍。

傳統中值濾波演算法的執行效率由快速排序演算法的時間複雜度決定，其基本思想是通過每次排序將像素點集合劃分為相互獨立的兩個子集合，其中一個子集合中所有像素值都比另一個子集合中所有像素值小，然後分別對這兩個子集合繼續進行排序，直到濾波窗口中所有像素點有序。

中值濾波是一種非線性平滑技術，它將每一像素點的灰度值設置為該點某鄰域窗口內的所有像素點灰度值的中值，是基於排序統計理論的一種能有效抑制雜訊的非線性信號處理技術，中值濾波的基本原理是把數字圖像或數字序列中一點的值用該點的一個鄰域中各點值的中值代替，讓周圍的像素值接近的真實值，從而消除孤立的雜訊點。方法是用某種結構的二維滑動模板，將板內像素按照像素值的大小進行排序，生成單調上升（或下降）的為二維數據序列。

在中值濾波窗口內各點有相同的輸出作用，若強調中間點或離該點較近的作用點，可改變窗口中變數個數，使多個變數值等於一個點值，再對擴展後的灰度值數字序列求中值。缺點是對邊緣像素進行窗口擴展後，將超出圖像邊界，引起邊界效應。偽代碼如下：

改進：僅存儲當前背景，如果當前幀在某個像素處變亮了，背景的亮度值就加1，如果變暗了，背景值就減1

基於混合高斯模型（GMM）的背景建模

高斯混合模型（Gaussian Mixed Model）指的是多個高斯分布函數的線性組合，理論上GMM可以擬合出任意類型的分布，通常用於解決同一集合下的數據包含多個不同的分布的情況（或者是同一類分布但參數不一樣，或者是不同類型的分布，比如正態分布和伯努利分布）。

基於混合高斯模型的背景建模將圖像中的每個像素看成是從混合高斯分布樣本中採樣得到的隨機變數，在這些高斯分布中，一些高斯分布表示背景，另一些表示前景（運動物體），根據演算法估計出每個像素點屬於哪一個高斯模型，進而判斷該像素是前景或背景；K的選擇由計算效率等因素決定，通常情況下，K的取值範圍為3~7；設特定像素在時刻t的亮度值為gt，用K個高斯分布表示像素gt的特徵，觀測概率

$p(g_t)=sum_{k=1}^K omega_{kt} frac{1}{sqrt{2pi}} exp left(-frac{(g_t-mu_{kt})^2}{sigma_{kt}^2} ight)\ sum_{k=1}^K omega_{kt}=1$

下面具體介紹演算法步驟。

第一步：模型初始化：確定高斯模型個數K，初始化模型各個參數

- [均值] 第一個高斯模型的均值取輸入視頻第一幀對應的像素值，其他設為0

- [方差] 初始方差相等，取值與該視頻的動態特性有關

- [權重ω] 在初始化時，一般第一個高斯模型的權重取較大值，其他取較小值

第二步：背景描述：確定描述背景模型的高斯函數；假設所有幀中，背景像素的比例超過閾值T，所有K個高斯函數按照表達式 $omega_{kt}/sigma_{kt}$ 排序，比值越大，該高斯函數代表背景的可能性越大，排序後，前B個高斯模型k=1,2,…,B作為背景 $B= ext{arg} min_b sum_{k=1}^b omega_{kt}>T$

第三步：前景判斷：對於任意像素，比較是否滿足各個高斯模型

$l_{kt}=egin{cases} 1 & |g_t-mu_{k(t-1)}| le lambda sigma_{k(t-1)} \ 0 & ext{others} end{cases}$

計算該像素滿足高斯模型的最小編號，判斷該像素是背景或前景。

第四步：模型更新：如果沒有找到匹配的高斯模型，將第K個(最後一個)高斯模型用一個新的高斯函數替換，新的高斯函數具有像素點的均值、一個較大的方差和較小的權重

$egin{cases} mu_{kt}=g_t \ sigma_{kt}^2 = 2max_k sigma_{k(t-1)}^2 \ omega_{kt} =0.5 min_k omega_{k(t-1)} end{cases} ag{1}$

重新歸一化各高斯模型的權重，或者採用梯度下降策略，獲取第t幀圖像，亮度向量為$x_t$，決定哪些高斯函數與之匹配，從中選擇最好的，記為l

$omega_{kt}= egin{cases} (1-eta)omega_{k(t-1)} & k e l \omega_{k(t-1)} & k = l end{cases} ag{2}$

$egin{cases} ho=eta N(x_t|mu_{l(t-1)},sigma_{l(t-1)}^2)\ mu_{lt}=(1- ho)mu_{l(t-1)}+ ho x_t\ sigma_{l(t-1)}^2=(1- ho)sigma_{l(t-1)}^2 + ho(x_t-mu_{l(t-1)})^T (x_t-mu_{l(t-1)}) end{cases}$

總結如下：

演算法：基於混合高斯模型的背景減除

初始化：選擇高斯函數數目和學習率η

repeat

獲取第t幀圖像，亮度向量為xt，決定哪些高斯函數與之匹配，從中選擇最好的，記為l，如果找到一個匹配高斯函數l，按式(2)更新所有高斯函數的權重和高斯函數l的參數

如果沒有高斯函數匹配像素，丟棄權重最低的高斯函數，用一個按照式(1)構造的新的高斯函數代替

根據高斯函數的權重，確定描述背景模型的高斯函數，判斷當前像素是背景還是前景

處理所有像素點，對前景圖像使用模糊以及形態學上的膨脹腐蝕組合操作去除小區域，填充大目標中的空洞，獲得場景中的運動目標

until 所有幀都被處理

貝葉斯濾波和粒子濾波

在卡爾曼濾波中，我們建立的模型

$egin{cases} X_t=A_{t,t-1}X_{t-1}+W_t\ Z_t=C_tX_t+V_t end{cases}$

具有這樣兩條性質

- [狀態量的馬爾可夫性] k時刻狀態量只與k-1時刻的狀態量有關

- [觀測量的條件獨立性] k時刻觀測值只與k時刻狀態值有關

如果我們把位置和觀測建模為概率分布，即

$egin{cases} x_k sim p(x_k|x_{1:k-1},y_{1:k-1})\ y_k sim p(y_k|x_{1:k},y_{1:k-1}) end{cases}$

根據貝葉斯定理

$p(x_{k}|y_{1:k})=frac{p(y_{1:k}|x_{k})p(x_{k})}{p(y_{1:k})}$

這種方法也就是所謂的批處理貝葉斯方法。但是這種方法每次都要拿所有的測量值來重新計算概率分布，對於計算機是個沉重的負擔。因此，有了改進的基於遞歸的貝葉斯濾波方法

$p(x_k|y_{1:k-1})=int{p(x_k|x_{k-1})p(x_{k-1}|y_{1:k-1})} dx_{k-1}$

其中 $p(x_k|x_{k-1})$ 代表了狀態由k-1時刻到k時刻的轉移概率。利用觀測值進行更新

$p(x_k|y_{1:k})=frac{p(y_k|x_k)p(x_k|y_{1:k-1})}{int{p(y_k|x_k)p(x_k|y_{1:k-1})} dx_k}$

這個後驗分布，將作為估計k+1時刻位置的先驗分布，從而開始下一輪的遞歸解算。

粒子濾波的結構實際上就是加入蒙特卡羅方法（Monte Carlo method，即以某時間出現的頻率來指代該事件的概率）的卡爾曼濾波，該方法的基本思想是用一組樣本（或稱粒子）來近似表示系統的後驗概率分布，然後使用這一近似的表示來估計非線性系統的狀態。採用此思想，在濾波過程中粒子濾波可以處理任意形式的概率，而不像卡爾曼濾波只能處理線性高斯分布的概率問題。粒子濾波的一大優勢也在於此。

粒子濾波的三個重要步驟為

- [粒子採樣] 從建議分布中抽取一組粒子

- [粒子加權] 根據觀測概率分布，重要性分布以及貝葉斯公式計算每個粒子的權值

- [估計輸出] 估計輸出系統狀態的均值協方差等

此外，為了應對粒子退化現象，還採用了重採樣等策略。

meanshift

本小節主要參考了 https://blog.csdn.net/dcrmg/article/details/52694557

Mean Shift(均值漂移)演算法是無參密度估計理論的一種，無參密度估計不需要事先知道對象的任何先驗知識，完全依靠訓練數據進行估計，並且可以用於任意形狀的密度估計，在某一連續點處的密度函數值可由該點鄰域中的若干樣本點估計得出。

Mean shift將特徵空間視為先驗概率密度函數，那麼輸入就被視為是一組滿足某種概率分布的樣本點，這樣一來，特徵空間中數據最密集的地方，對應於概率密度最大的地方，且概率密度的質心就可以被視為是概率密度函數的局部最優值，也就是要求的聚類中心。對於每一個樣本點，計算以它為中心的某個範圍內所有樣本點的均值，作為新的中心（這就是shift，即圓心的漂移），移動直至收斂。這樣每一輪迭代，中心都會向數據更密集的地方移動，直到最後穩定收斂到樣本的「質心」。

可以直觀理解為：在樣本空間中，任選一個點，然後以這個點為圓心，劃定一個圓形的區域。在此區域內的所有點以圓心為起點，產生N個向量，然後把這些向量都相加，再以向量的終點為圓心，劃定同樣半徑的圓形區域，執行同樣操作，如此迭代，直到收斂，如下圖所示。

給定d維空間Rd的n個樣本點 ,i=1,…,n,在空間中任選一點x，那麼Mean Shift（均值漂移）向量的基本形式定義為: $M_h=frac{1}{K} sum_{x_i in S_K} (x_i-x)$

Sk是一個半徑為h的高維球區域，即滿足以下關係的y點的集合

$S_h (x)={ y:(y-x_i )^T (y-x_i )}$

其中k表示在這n個樣本點xi里，有k個點落入Sk區域中；考慮軟間隔優化問題：

$egin{aligned} min_{gamma,w,b} & , frac{1}{2} ||oldsymbol{w}||^2 +Csum_{i=1}^m xi_i\ ext{s.t.} & , y^{(i)} [w^T x^{(i) }+b] ge 1-xi_i quad xi_i ge 0 end{aligned}$

引入非負參數ξ後（稱為鬆弛變數），就允許某些樣本點的函數間隔小於1，即在最大間隔區間裡面，或者函數間隔是負數，即樣本點在對方的區域中。而放鬆限制條件後，我們需要重新調整目標函數，以對離群點進行處罰，這就是正則化方法。離群點越多，目標函數值越大，而我們要求的是儘可能小的目標函數值。這裡的C是離群點的權重，C越大表明離群點對目標函數影響越大，也就是越不希望看到離群點。我們看到，目標函數控制了離群點的數目和程度，使大部分樣本點仍然遵守限制條件。拉格朗日運算元：

$L(w,b,xi ,alpha ,r)=frac{1}{2} w^T w+Csum_{i=1}^m xi _i -sum_{i=1}^m alpha _i [y^{(i) } (w^T x^{(i)}+b)-1+xi _i ] -sum_{i=1}^m r_i xi _i$

這裡的α和γ都是拉格朗日乘子，結果如下：

$egin{aligned} max_{alpha} & , ?W(alpha )=sum_{i=1}^m alpha _i -frac{1}{2} sum_{i,j=1}^m y^{(i) } y^{(j) } alpha _i alpha _j<x^{(i) },x^{(j)}>\ ext{s.t.} & , alpha_i in [0,C] quad sum_{i=1}^m alpha _i y^{(i) }=0 end{aligned}$

對於含有不等式約束的優化問題，需要考慮KKT條件：

$egin{cases} alpha _i=0\ alpha _i=C\ alpha _iin (0,C) end{cases}Rightarrow y^{(i) } [w^T x^{(i) }+b]egin{cases} geq 1 \ leq 1 \ =1 end{cases}$

把基本的meanshift向量加入核函數k(x)，則

$f_{h,k} (x)=frac{C_{k,d}}{nh^d } sum_{i=1}^n frac{k(x)}{h^2} ||x-x_i||^2$

其中h為半徑， $C_{k,d}/n_{hd}$ 為單位密度，要使得上式f得到最大，對上式進行求導，梯度

$abla f=frac{2c_{k,d}}{nh^{d+2}} sum_{i=1}^n frac{(x-x_i )k(x)}{h^2} ||x-x_i||^2 =0$

則

$x=frac{sum_{i=1}^n x_i k ||x-x_i||^2}{sum_{i=1}^n k ||x-x_i||^2}$

下面介紹meanshift演算法怎樣運用到圖像上的聚類核跟蹤。一般一個圖像就是個矩陣，像素點均勻的分布在圖像上，就沒有點的稠密性。所以怎樣來定義點的概率密度，這才是最關鍵的。如果我們計算點x的概率密度，採用的方法如下：以x為圓心，以h為半徑。落在球內的點位xi ，認為

- x像素點的顏色與xi像素點顏色越相近，我們定義概率密度越高。

- 離x的位置越近的像素點xi，定義概率密度越高。

所以定義總的概率密度，是兩個規則概率密度乘積的結果，可以表示為

$K_{s,r} (x)=frac{C}{h_s^2 h_r^2 } k left( frac{1}{h_s^2} ||(x^{(s) }-x_i^{(s) }||^2 ight)cdot k left( frac{1}{h_r^2} ||(x^{(r) }-x_i^{(r) }||^2 ight)$

將meanshift演算法擴展到連續圖像序列，就是camshift演算法。它將視頻的所有幀做meanshift運算，並將上一幀的結果，即搜索窗的大小和中心，作為下一幀meanshift演算法搜索窗的初始值；如此迭代下去，就可以實現對目標的跟蹤。