馬 龍，王魯平，李 飚，陳小天

（國(guó)防科技大學(xué)ATR重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙410073）

0 引 言

光流技術(shù)在運(yùn)動(dòng)檢測(cè)、三維重建、機(jī)器人導(dǎo)視等領(lǐng)域有著重要應(yīng)用［1－3］，微分光流法是最重要的光流方法之一［4－8］。微分光流法是由 Horn和Schunck最先提出的［4］，他們?cè)诠饬骰痉匠痰幕A(chǔ)上引入了全局平滑約束，從而使光流問(wèn)題可適化，進(jìn)而推導(dǎo)出了光流估計(jì)式。微分光流法計(jì)算簡(jiǎn)單、幾何意義直觀，但有較大的局限性，突出表現(xiàn)在兩個(gè)方面：一是光流基本方程是基于灰度恒定假設(shè)建立的，在光照強(qiáng)度變化的環(huán)境下具有很大誤差；二是全局平滑約束會(huì)使光流場(chǎng)平滑越過(guò)物體邊界，造成光流場(chǎng)過(guò)平滑。

為提高光流方法的光變魯棒性，人們?cè)趥鹘y(tǒng)基本方程的基礎(chǔ)上推導(dǎo)了可適應(yīng)光變環(huán)境的基本方程，提出了許多新的方法［9－10］。其中，文獻(xiàn)［9］的方法最具代表性，文中假設(shè)光強(qiáng)發(fā)生線性變化并據(jù)此修正傳統(tǒng)方程，從而得到了新的基本方程?？稍诠庾儹h(huán)境下獲得較好的光流場(chǎng)估計(jì)結(jié)果。

為解決光流估計(jì)的過(guò)平滑問(wèn)題，人們提出了許多新的光流場(chǎng)平滑方法，比較有代表性的有Lucas-Kanade局部平滑法［5］、Nagel有向平滑法［6］和Black-Anandan分段平滑法［7］。其中，采用了最小二乘（Least squares，LS）技術(shù)的局部平滑法效果最優(yōu)［11］。但是，局部平滑方法也具有一定的局限性，突出表現(xiàn)為在物體邊緣處該方法的光流估計(jì)誤差很大。針對(duì)這一問(wèn)題，一些學(xué)者提出在局部平滑約束的框架下采用最小二乘中值方法（Least median of squares，LMedS）代替LS法求解光流［10］。這一改進(jìn)使得局部光流法在物體邊緣的光流估計(jì)精度得到很大提高，但是，在物體邊角位置，用LMedS法仍然不能得到準(zhǔn)確的光流估計(jì)。同時(shí)，LMedS法對(duì)噪聲較敏感。

本文首先從理論上對(duì)文獻(xiàn)［9］的新基本方程進(jìn)行了推導(dǎo)并討論了相關(guān)參量的物理含義，而后借助u-v平面解析了傳統(tǒng)局部方法的局限性，最后提出了一種新的基于像素篩選技術(shù)的光流估計(jì)方法，該方法不僅具有光變和噪聲魯棒性，而且能在物體邊角位置獲得準(zhǔn)確的光流估計(jì)。

1 光變環(huán)境下光流基本方程的推導(dǎo)

如圖1所示，設(shè)P為一物體表面元在t時(shí)刻的位置，t＋δt（δ符號(hào)表示變量的微小變化）時(shí)刻該表面元運(yùn)動(dòng)到位置P′。在此期間，表面元在像平面上的投影由像素p移動(dòng)到像素p′，p和p′的像平面坐標(biāo)分別設(shè)為（x，y）和（x＋δx，y＋δy）。又設(shè)t時(shí)刻光照強(qiáng)度為L(zhǎng)（t）＝l（t）1L，其中l(wèi)（t）為t時(shí)刻光強(qiáng)的大小，1L為光強(qiáng)單位方向矢量。則像素p的灰度值

假設(shè)t＋δt時(shí)刻光強(qiáng)發(fā)生了變化，不失一般性，設(shè)此變化是線性的，即l（t＋δt）＝（1＋δk）l（t）＋δc。其中，δk和δc分別為光強(qiáng)的乘性和加性增量，它們?yōu)闀r(shí)間t的函數(shù)。則p′的灰度值

圖1 運(yùn)動(dòng)物體成像灰度變化示意圖Fig.1 Sketch map of intensity variance for moving object

式（1）（2）是假設(shè)物體表面為朗伯面而得到的，其中，ρ為表面元漫反射系數(shù)，為一常量；1P和1P′分別為表面元在位置P和P′時(shí)的單位法向量，在δt內(nèi)可認(rèn)為1P＝1P′。

由式（1）（2）可知，在δt時(shí)間內(nèi)，表面元成像灰度變化量為

上式中的第一項(xiàng)與光強(qiáng)相關(guān)，其中ρl （t） 1P·1L＝I （x ，y，t）；第二項(xiàng)僅與光強(qiáng)方向有關(guān)，而與光強(qiáng)大小無(wú)關(guān)。由于假設(shè)光強(qiáng)只發(fā)生線性變化，即1L為時(shí)間不變量，因此ρ1P·1L也為時(shí)間不變量，設(shè)為m1。

式（3）左端的I （x＋δx，y＋δy，t＋δt） 根據(jù)Taylor公式展開(kāi)得

式中：Ix（ x ，y，t）、Iy（ x ，y，t）和It（ x ，y，t）分別表示I （x ，y，t）的位置和時(shí)間導(dǎo)數(shù)；O是關(guān)于δx、δy、δt的二階及二階以上的項(xiàng)。將式（4）代入式（3）左端并化簡(jiǎn)得

式（5）兩端同除以δt，而后取極限δt→0，得：

又令m′＝m1m2，則

式（8）即為光變環(huán)境下的光流基本方程。這一方程是在考慮光強(qiáng)變化的情況下得出的，在光強(qiáng)恒定的情況下，光強(qiáng)乘性變化率m和加性變化率m2均為0，此時(shí)，該方程退化為傳統(tǒng)基本方程。因此，式（8）可以看作是傳統(tǒng)基本方程在光變環(huán)境下的推廣。

2 經(jīng)典局部光流法的性能分析

局部光流法求解光流場(chǎng)的一般思路為：假定光流場(chǎng)是局部平滑的，即當(dāng)前像素x的鄰域Ω內(nèi)（一般為N×N鄰域窗口，設(shè)n＝N×N）所有像素的光流值相同，由此可得超定方程組

求解此超定方程組即可得到x處的光流估計(jì)。

記殘差ri（u，v）＝｜Ix（xi，yi，t）u＋I(xiàn)y（xi，yi， t）v＋I(xiàn)t（xi，yi，t）－mI（xi，yi，t）－m′｜，則用LS法求解光流的估計(jì)式為

LMedS法求解光流的估計(jì)式為

超定方程組（9）可看作是u-v平面上的一簇直線，求解超定方程組即是估計(jì)這簇直線在u－v平面上的交點(diǎn)位置。

當(dāng)Ω落于物體內(nèi)部區(qū)域時(shí)，Ω中所有像素的光流值近似相等（見(jiàn)圖2（a），三角形標(biāo)記的為當(dāng)前像素），直線簇近似相交于u-v平面上一點(diǎn)（u，v），則（u，v）即為超定方程組的解，也即像素x的光流值（見(jiàn)圖2（b）中“十”字標(biāo)記），此時(shí)，LS法和LMedS法均可求得光流的精確估計(jì)。

圖2 傳統(tǒng)局部光流方法的u-v平面解析Fig.2 u-v analyses for traditional local flow method

當(dāng)Ω落在物體邊緣時(shí)，Ω中可能同時(shí)包含物體內(nèi)部像素和外部像素。設(shè)內(nèi)部像素和外部像素的光流分別為 （u′，v′）和 （u″ ，v″）（圖2中分別用“十”字和“叉”號(hào)標(biāo)記），且 （u′，v′）≠（u″，v″）。此時(shí)，對(duì)應(yīng)內(nèi)部像素和外部像素的兩組直線簇在u-v平面上有兩個(gè)不同的公共交點(diǎn)。在這種情況下，LS法估計(jì)的光流值在圖2（d）中用實(shí)三角標(biāo)記，顯然這一估計(jì)值與實(shí)際光流值 （u′ ，v′）的偏差很大。對(duì)于LMedS法而言，如果Ω中內(nèi)部像素?cái)?shù)量占優(yōu)（見(jiàn)圖2（c）），則可根據(jù)式（11）獲得準(zhǔn)確的光流估計(jì)（見(jiàn)圖2（e）實(shí)三角標(biāo)記）；否則，根據(jù)式（11）將不能獲得準(zhǔn)確的光流估計(jì)（見(jiàn)圖2（f）實(shí)三角標(biāo)記）。

由以上分析可知，在物體內(nèi)部區(qū)域，用LS法和LMedS法都可以得到準(zhǔn)確的光流估計(jì)；在物體邊緣，當(dāng)Ω中的內(nèi)部像素?cái)?shù)量不占優(yōu)時(shí)，兩種方法都不能獲得正確的光流估計(jì)。而在實(shí)際光流估計(jì)過(guò)程中，當(dāng)鄰域窗口落于物體邊角位置時(shí)，Ω中內(nèi)部像素少于外部像素的情況會(huì)經(jīng)常出現(xiàn)。因此，有必要探索能對(duì)物體邊緣處光流進(jìn)行正確估計(jì)的新方法。

3 本文方法

LMedS法只適用于Ω中物體內(nèi)部像素?cái)?shù)量占優(yōu)時(shí)的光流場(chǎng)估計(jì)，當(dāng)Ω落于物體邊角位置時(shí)，可以設(shè)法舍棄一部分外部像素，從而使LMedS法的適用條件得以滿足。這需要同時(shí)考慮兩個(gè)問(wèn)題：①Ω中所有像素的歸屬都是未知的，因而需要確定哪些像素是可以舍棄的外部像素；②Ω中內(nèi)部像素的數(shù)量是未知的，因而需要確定應(yīng)該舍棄掉多少個(gè)外部像素。

問(wèn)題①是一個(gè)模式分類問(wèn)題，由于問(wèn)題②并不要求舍棄所有的外部像素，因此，不需要對(duì)圖像進(jìn)行精確分割，只需要把Ω內(nèi)最可能歸于外部區(qū)域的像素標(biāo)記出來(lái)即可。

對(duì)于問(wèn)題②，假設(shè)物體成像表面是凸的，且不存在線狀部分（成像表面不存在兩個(gè)以上連續(xù)的單像素點(diǎn)）。則當(dāng)鄰域窗口Ω落于該物體上時(shí)，Ω內(nèi)物體內(nèi)部像素至少為2 N－1個(gè)。因此，應(yīng)舍棄的像素?cái)?shù)k＝N2－（2 N－1）－（2 N－2）＝N2－4 N－3個(gè)。事實(shí)上，通常情況下，Ω中物體內(nèi)部像素多于2 N－1個(gè)，因此被舍棄的k個(gè)像素中可能包含內(nèi)部像素，但這不會(huì)影響Ω中物體內(nèi)部像素占優(yōu)的處理結(jié)果。

基于以上分析，提出了一種基于像素篩選技術(shù)的光流估計(jì)方法（Pixels-filtrated-based method，PFBM）。這種方法首先對(duì)Ω中的像素進(jìn)行篩選，使得物體內(nèi)部像素占優(yōu)，并在此基礎(chǔ)上采用LMedS估計(jì)初始光流；而后計(jì)算初始光流在窗口內(nèi)每一像素處的殘差值，并根據(jù)殘差值大小對(duì)窗口內(nèi)像素進(jìn)行二次篩選；最后，以篩選后的保留像素為基礎(chǔ)，采用最小二乘法估計(jì)最終光流值。

3.1 像素初次篩選及初始光流估計(jì)

像素初次篩選的目的是舍棄一定數(shù)量的外部像素，使得Ω中內(nèi)部像素?cái)?shù)量占優(yōu)。這里的內(nèi)部像素是指與當(dāng)前像素屬于同一物體的像素。對(duì)于灰度圖像而言，判斷兩個(gè)像素是否屬于同一成像物體的較為方便、可靠的方法是度量其像素灰度的接近程度?；谶@一分析，建立中心灰度差分矩陣來(lái)度量Ω中各像素屬于內(nèi)部像素的可能程度。

式中：Ta為一自適應(yīng)閾值，它使得MΩ中0元素?cái)?shù)為N2－4 N＋3。

經(jīng)過(guò)篩選，Ω中的內(nèi)部像素?cái)?shù)量占優(yōu)，此時(shí)，由LMedS法可得

上式是LMedS的加權(quán)形式，稱之為最小加權(quán)二乘中值（Least weighted median of squares，LWMS）方法。

3.2 像素二次篩選及最終光流估計(jì)

在Ω中像素初次篩選的基礎(chǔ)上，采用LWMS法得到了初始光流估計(jì)（u′，v′）。但同LMedS法一樣，LWMS法的光流估計(jì)也是建立在單個(gè)殘差值基礎(chǔ)上的，因此不具有噪聲魯棒性。為提高算法的抗噪性能，把（u′，v′）作為初始光流，并在此基礎(chǔ)上估計(jì)最終光流值。

由3.1節(jié)可知，（u′，v′）為Ω中內(nèi)部像素的共同光流值，因此，內(nèi)部像素在（u′，v′）上的殘差值應(yīng)該較小，而外部像素在（u′，v′）上的殘差值應(yīng)該較大。據(jù)此，可對(duì)Ω中像素進(jìn)行二次篩選，把所有可能的內(nèi)部像素挑選出來(lái)。具體采用Rousseeuw等［12］提出的方法，即

（1）根據(jù)ri（u′，v′）＝｜Ix（xi，yi，t）u′＋I(xiàn)y（xi，yi，t）v′＋I(xiàn)t（xi，yi，t）－mI（xi，yi，t）－m′｜（i＝0，…，n－1）計(jì)算（u′，v′）在Ω上的殘差，并計(jì)算初始?xì)埐铋撝?/p>

式中：n為Ω中的像素個(gè)數(shù)；p為所需估計(jì)的參量個(gè)數(shù)，這里取為4。

（2）對(duì)Ω上的像素進(jìn)行篩選，得到初始加權(quán)系數(shù)

并計(jì)算殘差閾值σ

（3）用σ替代σ′代入式（15），對(duì)Ω上的像素進(jìn)行重新篩選，得到加權(quán)系數(shù)wi（i＝0，…，n－1）。

二次篩選后的保留像素被認(rèn)為是物體內(nèi)部像素，最終光流值以保留像素為基礎(chǔ)，采用最小加權(quán)二乘方法得到

4 實(shí)驗(yàn)與分析

采用邊緣較為明顯的Marbled block圖像序列進(jìn)行算法性能測(cè)試。光流估計(jì)精度通過(guò)角度誤差來(lái)度量，角度誤差定義為

式中：（uc，vc）表示實(shí)際光流場(chǎng)；（ue，ve）表示估計(jì)光流場(chǎng)；ψ的均值和標(biāo)準(zhǔn)差分別用ψ－和σψ表示。

4.1 物體邊緣處的光流估計(jì)實(shí)驗(yàn)

在圖3所示的Marbled block圖像中，分別在大理石塊和離其最近的立方柱上各選擇一個(gè)尖角位置和邊緣位置，在所選4個(gè)像素的鄰域窗口內(nèi)Ω（5×5），物體內(nèi)部像素?cái)?shù)目均不占優(yōu)（見(jiàn)圖3）。表1為PFBM方法的估計(jì)誤差，同時(shí)給出了對(duì)比算法的估計(jì)誤差。

圖3 Marbled block圖像Fig.3 Marbled block image

從表1可以看出，6種光流方法中，PFBM方法無(wú)論是在大理石塊還是在立方柱的尖角和邊緣位置都得到了較好的光流估計(jì)。

表2為6種光流方法對(duì)原始圖像的光流場(chǎng)估計(jì)誤差，總體而言，PFBM法的估計(jì)誤差最小，這是因?yàn)镻FBM法在物體尖角和邊緣處得到的光流估計(jì)較為準(zhǔn)確，從而使得總體光流估計(jì)性能最優(yōu)。

表1 尖角邊緣處光流估計(jì)的角度誤差Table 1 Angle errors of estimated flow in corners and edges

表2 光流場(chǎng)估計(jì)誤差Table 2 Errors of estimated flow field

4.2 抗光變及抗噪聲實(shí)驗(yàn)

為測(cè)試算法的光變魯棒性，首先采用高斯光變模型對(duì)Marbled block圖像進(jìn)行灰度變換以模擬光照條件的變化對(duì)圖像灰度的影響，而后分別采用6種光流方法估計(jì)光照變化后的圖像光流場(chǎng)，所得結(jié)果如表2所示。從實(shí)驗(yàn)結(jié)果可以看出，K-LMedS法和PFBM法的光流估計(jì)誤差較小，而其他方法的估計(jì)誤差則相當(dāng)大。這是因?yàn)镵LMedS法和PFBM法考慮了光照變化，并且在其光流基本方程中加入了光變約束。

為測(cè)試算法的抗噪性能，對(duì)Marbled block圖像分別加入了4個(gè)等級(jí)的高斯噪聲，各種方法對(duì)噪聲圖像的光流估計(jì)誤差如表2所示。實(shí)驗(yàn)結(jié)果中，LMedS、K-LMedS、PFBM方法給出的光流估計(jì)誤差較小，這是因?yàn)檫@幾種方法最終的光流求解都是基于最小二乘法進(jìn)行的。噪聲對(duì)光流場(chǎng)估計(jì)誤差的影響如圖4所示。

圖4 噪聲對(duì)光流估計(jì)性能的影響Fig.4 Noise impact on performance of flow estimation

5 結(jié)束語(yǔ)

提出了一種基于像素篩選技術(shù)的光流估計(jì)方法，該方法可以有效解決物體邊緣的光流估計(jì)問(wèn)題。由于采用了光變條件下的基本方程和LS技術(shù)，該算法對(duì)光照變化及噪聲干擾具有較強(qiáng)的魯棒性。

［1］Axel Besinger，Tamara Sztynda，Saroj Lal，et al. Optical flow based analyses to detect emotion from human facial image data［J］.Expert Systems with Applications，2010，37（12）：8897-8902.

［2］Lee Dah-Jye，Merrell Paul，Wei Zhao-yi，et al. Two-frame structure from motion using optical flow probability distributions for unmanned air vehicle obstacle avoidance［J］.Machine Vision and Applications，2010，21（3）：229-240.

［3］Kai Guo，Ishwar P，Konrad J，et al.Action recognition using sparse representation on covariance manifolds of optical flow［C］∥IEEE International Conference on AVSS，Boston，MA，2010：188-195.

［4］Horn B K，Schunck B，Schunck Brian G.Determining optical flow［J］.Artificial Intelligence，1981，17（1）：185-203.

［5］Lucas B D，Kanade T.An iterative image registration technique with an application to stereo vision［C］∥Proc DARPA Image Understanding Workshop，1981：121-130.

［6］Nagel H H，Enkelmann W.An investigation of smoothness constraints for the estimation of displacement vector field from image sequences［J］. IEEE Transactions on Pattern Analyze and Machine Intelligence，1986，8：565-593.

［7］Black M J，Anandan P.The robust estimation of multiple motions：parametric and piecewise smooth flow fields［J］.Computer Vision and Image Understanding，1996，63：75-104.

［8］Brox Thomas，Malik Jitendra.Large displacement optical flow：descriptor matching in variational motion estimation［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2011，33（3）：500-513.

［9］Kim Yeon-Ho，Kak Avinash C.Error analysis of robust optical flow estimation by least median of squares method for the varying illumination model［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2006，28（9）：1418-1435.

［10］Lai Shang-Hong，F(xiàn)ang Ming.Robust and efficient image alignment with spatially varying illumination models［C］∥Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition，1999：167-172.

［11］Barron J L，F(xiàn)leet D J，Beauchemin S S.Performance of optical flow techniques［J］.International Journal of Computer Vision，1994，12（1）：43-77.

［12］Rousseeuw Peter J.Least median of squares regression［J］.Journal of the American Statistical Association，1984，79（388）：871-880.