鄭 佳， 王洪雁， 裴炳南

(大連大學(xué)遼寧省北斗高精度位置服務(wù)技術(shù)工程實(shí)驗(yàn)室，遼寧 大連 116622)

0 引言

流體力學(xué)及空氣動(dòng)力學(xué)研究中流體速度場(chǎng)測(cè)量對(duì)于了解復(fù)雜流體具有重大意義。流體運(yùn)動(dòng)是一種典型的非剛性運(yùn)動(dòng)，其計(jì)算須基于圖像處理技術(shù)?？赏ㄟ^(guò)對(duì)運(yùn)動(dòng)圖像序列分析獲得局部流體運(yùn)動(dòng)矢量大小、方向及分布情況，進(jìn)而獲取諸如粘性及渦流場(chǎng)分布等物理特性[1]。通常情況下，由于運(yùn)動(dòng)物體透明或不易通過(guò)光學(xué)設(shè)備觀測(cè)，需將可見(jiàn)粒子置入被測(cè)物，以通過(guò)估計(jì)粒子運(yùn)動(dòng)矢量間接獲得流體運(yùn)動(dòng)特征，此即謂粒子圖像測(cè)速[2](Particle Image Velocimetry，PIV)。

常見(jiàn)光流計(jì)算方法，如變分光流方法的分辨率可達(dá)到亞像素級(jí)，同時(shí)，假設(shè)流體連續(xù)時(shí)空變化特性與光流方程中圖像序列局部時(shí)空可微本質(zhì)相同，因而，基于光流的PIV得到了廣泛的關(guān)注。

馬鵬飛等提出基于Lucsa-Kanade(LK)局部光流的PIV估計(jì)方法[3]，此算法計(jì)算復(fù)雜度低，穩(wěn)健性較好，然而光流場(chǎng)邊界較為模糊，邊緣像素點(diǎn)光流估計(jì)較差，從而所得光流場(chǎng)較為稀疏；針對(duì)此問(wèn)題，孫立志等提出基于金字塔LK的多尺度光流算法以計(jì)算PIV[4]，所提算法可解決在運(yùn)動(dòng)目標(biāo)速度過(guò)大、相鄰幀連續(xù)性不強(qiáng)時(shí)的光流計(jì)算問(wèn)題，所得光流精度高于LK算法，然而得到的光流場(chǎng)仍較為稀疏；為此，黃湛等提出基于Horn-Schunck(HS)的PIV計(jì)算方法[5]，此算法可得稠密光流場(chǎng)，然而邊界容易模糊，計(jì)算復(fù)雜度高，且穩(wěn)健性較差；為進(jìn)一步改善光流估計(jì)精度，余學(xué)敏等使用全局與局部結(jié)合的方法，獲得了更高的光流估計(jì)精度[6]，然而此方法所得光流在邊緣地區(qū)穩(wěn)健性較差；為改善光流估計(jì)穩(wěn)健性，文獻(xiàn)[7]在局部與全局結(jié)合的基礎(chǔ)上加入各向異性濾波器以減少邊緣光流擴(kuò)散，但是此方法對(duì)噪聲及異常點(diǎn)的穩(wěn)健性較差；文獻(xiàn)[8]在全局與局部結(jié)合的基礎(chǔ)上加入懲罰項(xiàng)，以增強(qiáng)對(duì)噪聲和異常點(diǎn)的穩(wěn)健性；文獻(xiàn)[9]基于物理學(xué)的光流方法從可視化圖像獲取高分辨速度場(chǎng)，然而由于邊緣區(qū)域光流擴(kuò)散及噪聲和異常點(diǎn)的影響，光流穩(wěn)健性較差。針對(duì)上述問(wèn)題，本文提出一種基于物理學(xué)的穩(wěn)健光流計(jì)算方法以改善光流計(jì)算穩(wěn)健性。在基于物理學(xué)的光流方法中引入各向異性濾波器以增強(qiáng)邊緣光流的穩(wěn)健性，增加懲罰因子以減少噪聲及異常點(diǎn)對(duì)光流計(jì)算的影響，而后基于變分方法極小化光流能量函數(shù)以求解歐拉-拉格朗日方程，最后通過(guò)迭代方法求得速度場(chǎng)。

1 經(jīng)典光流法

1.1 光流約束方程

在相鄰圖像時(shí)間間隔和圖像灰度變化均很小的假設(shè)下，HORN等提出如下灰度光流計(jì)算思想：物體運(yùn)動(dòng)會(huì)引起與之對(duì)應(yīng)像素亮度連續(xù)變化，從而形成連續(xù)變化光流場(chǎng)[5]，其數(shù)學(xué)模型如下。

設(shè)圖像中一點(diǎn)(x，y)在t時(shí)刻的亮度為I(x，y，t)，在Δt時(shí)刻亮度為I(x+Δx,y+Δy,t+Δt)，當(dāng)Δt無(wú)窮小時(shí)可認(rèn)為亮度不變，則可得

I(x,y,t)=I(x+Δx,y+Δy,t+Δt)

(1)

對(duì)上式泰勒級(jí)數(shù)展開(kāi)，當(dāng)Δt→0可得

(2)

Ixu+Iyv+It=0

(3)

式中，u，v表示速度矢量的二維分量。由此可得光流計(jì)算基本等式。需要注意的是，式(3)中有兩個(gè)未知參數(shù)，因而此問(wèn)題無(wú)法求解。為獲取兩個(gè)速度場(chǎng)參數(shù)的估計(jì)，HORN等提出了基于全局約束的方法，而KANADE等則提出基于局部約束的方法。

1.2 Lucas-Kanade光流法

此方法核心思想如下：假設(shè)像素a與相鄰區(qū)域內(nèi)所有像素光流矢量相同，對(duì)區(qū)域內(nèi)像素賦予不同權(quán)重，則光流計(jì)算可等價(jià)為最小化如下能量函數(shù)[10]，即

(4)

式中：Ω表示以a點(diǎn)為中心的小區(qū)域；▽表示梯度算子；W(x)為窗函數(shù)，表示區(qū)域中像素mi(i=1，2，…，n)的權(quán)重，離a點(diǎn)越近，權(quán)重越高。上式可等價(jià)為如下最小二乘形式，即

A2W2Ap=ATW2b

(5)

式中，

(6)

diag(·)為對(duì)角陣。由式(5)可得

p=(ATW2A)-1ATW2b。

(7)

LK方法可得到稀疏光流場(chǎng)，相對(duì)于下述HS算法，計(jì)算復(fù)雜度低，然而穩(wěn)健性較差[3]。

1.3 Horn-Schunck光流法

HORN等提出的約束條件則是極小化平滑約束項(xiàng)[11]Es，即

(8)

而基于光流基本等式可得光流誤差Ec為

Ec=(Ixu+Iyv+It)2dxdy。

(9)

基于上述兩個(gè)約束條件可知，光流矢量可通過(guò)如下最小化問(wèn)題獲得，即

E={(Ixu+Iyv+It)2+λ(‖▽u‖2+‖▽v‖2)}dxdy

(10)

此問(wèn)題可通過(guò)基于梯度的方法進(jìn)行求解，從而降低計(jì)算復(fù)雜度，且可得到高精度像素瞬時(shí)位置速度及稠密光流，然而，此方法分辨率較低，邊緣光流易擴(kuò)散，且穩(wěn)健性較差[5]。

2 基于物理的光流法

針對(duì)上述算法分辨率較低的問(wèn)題，文獻(xiàn)[9]提出一種基于物理學(xué)的高分辨率光流計(jì)算方法。此算法詳細(xì)地給出了光流的物理意義，即光流與可視化圖像中流體或者顆粒的速度成正比[12]。此外，在圖像坐標(biāo)系中，物體空間坐標(biāo)系可通過(guò)透視投影變換來(lái)表述，所有這些情況下的投影運(yùn)動(dòng)方程可用基于物理的光流方程表述[9]為

(11)

當(dāng)g▽·p=0及f=0時(shí)，基于物理學(xué)的光流方程可以簡(jiǎn)化為基于HS光流算法的亮度約束方程[9]，即

(12)

極小化式(12)可得歐拉-拉格朗日方程[9]為

(13)

式(13)可通過(guò)標(biāo)準(zhǔn)有限差分方法來(lái)進(jìn)行光流求解?；谖锢韺W(xué)的光流計(jì)算方法，可以得到較HS方法更高分辨率的光流。然而，由于基于物理光流方法光流在邊緣地區(qū)易擴(kuò)散，且易受噪聲和異常點(diǎn)影響，光流整體穩(wěn)健性較差[9]。

3 基于物理學(xué)的穩(wěn)健光流算法

為了減少角落和邊緣光流擴(kuò)散對(duì)光流計(jì)算的影響，文獻(xiàn)[7]提出各向異性擴(kuò)散的光流計(jì)算方法，并給出了擴(kuò)散系數(shù)方程。此方法可減少光流在邊緣的傳播，增強(qiáng)邊緣光流穩(wěn)健性。所提的兩個(gè)擴(kuò)散系數(shù)方程可表示為

(14)

(15)

其中，式(14)更適用于高對(duì)比度情況下的穩(wěn)健光流求解，式(15)更適用于較寬區(qū)域的穩(wěn)健光流求解[13]。

常量K控制邊緣的敏感度，通常由實(shí)驗(yàn)確定其數(shù)值。為了能夠兼顧對(duì)比度和區(qū)域大小對(duì)光流計(jì)算的影響，求得穩(wěn)健性更高的光流，本文采用擴(kuò)散系數(shù)方程[7]為

D(‖▽I‖)=e-α‖▽I‖β1。

(16)

基于式(16)，光流能量函數(shù)式(12)可改寫(xiě)為

(17)

式中，E2=‖D▽u‖2+‖D▽v‖2。

加入懲罰項(xiàng)，式(17)可改寫(xiě)為

(18)

根據(jù)變分原理將E(u,v)最小化，可得歐拉-拉格朗日方程[14]為

(19)

式中：

L=L(u,v,ux,uy,vx,vy)=(E1)2+E3+λE2；

(20)

(21)

將式(20)及式(21)代入式(19)可化簡(jiǎn)為

(22)

式中，▽·(gp)=g▽·p+p·▽g，可用某點(diǎn)與其周圍速度平均值之差近似表示，即

(23)

結(jié)合式(23)，式(18)可重新表示為

(24)

式中，E4=gt+ugx+vgy+gux+gvy-f。

對(duì)式(24)整理可得

(25)

式中，

(26)

求解式(25)得un+1和vn+1為

(27)

式中，

(28)

至此，可得基于物理的穩(wěn)健光流計(jì)算方法。此方法可具體描述如下。

1) 讀取連續(xù)兩幀粒子圖像，對(duì)其濾波以減小噪聲，并初始化所需參數(shù)。

2) 使用sobel算子對(duì)濾波后圖像求x方向、y方向?qū)?shù)；利用兩幀中對(duì)應(yīng)像素相減方式求t方向?qū)?shù)。

3) 采用九點(diǎn)差分格式計(jì)算u，v均值[15]。

4) 計(jì)算求解點(diǎn)灰度梯度，設(shè)定速度平滑權(quán)重系數(shù)(設(shè)為1)，初始速度設(shè)為0。

5) 根據(jù)式(26)，基于Gauss-Seidel方法迭代求解[16]。

6) 若相鄰兩次迭代光流之差小于給定閾值，或迭代次數(shù)超過(guò)給定值，則終止迭代。

4 實(shí)驗(yàn)仿真及分析

通過(guò)與經(jīng)典HS，LK光流算法以及金字塔LK光流算法和基于物理的光流算法進(jìn)行對(duì)比，驗(yàn)證所提基于物理的穩(wěn)健算法針對(duì)所得光流密度及邊緣光流擴(kuò)散效果的有效性。仿真參數(shù)設(shè)置如下：各向異性濾波器相關(guān)參數(shù)α=5，β1=0.8，懲罰函數(shù)相關(guān)參數(shù)β2=0.02，迭代誤差ε=1×10-3，迭代次數(shù)n=200。

(29)

式中，

(30)

實(shí)驗(yàn)1 圖1所示為等流速狀態(tài)下兩幀連續(xù)的粒子圖像，以下算法皆基于此計(jì)算光流。圖2所示為L(zhǎng)K光流算法、HS光流算法、基于金字塔LK光流算法以及基于物理算法得到的光流圖像。由圖2可知，與HS光流算法相比，LK光流算法得到的光流稀疏，穩(wěn)健性較好。HS光流算法雖可得稠密光流，但穩(wěn)健性差。較HS光流算法和LK光流算法，基于金字塔LK光流算法所得光流精度較高，穩(wěn)健性較強(qiáng)，但光流較稀疏。此3種光流算法所得光流分辨率較低?；谖锢砉饬魉惴傻酶叻直婀饬?，但邊緣易擴(kuò)散，且易受噪聲和異常點(diǎn)影響，光流較稀疏，穩(wěn)健性較差。

圖1 原始圖像Fig.1 Original images

圖2 基于物理光流算法與經(jīng)典光流算法所得光流對(duì)比

實(shí)驗(yàn)2 圖3分別為基于物理光流算法及加入懲罰因子后基于物理光流算法所得光流圖像，懲罰因子中β2=0.02。由圖3可知，在基于物理光流算法基礎(chǔ)上加入懲罰因子后對(duì)噪聲和異常點(diǎn)的穩(wěn)健性較好，但由于光流在邊緣地區(qū)易擴(kuò)散，所得光流較為稀疏。

圖3 懲罰因子對(duì)基于物理光流算法的影響Fig.3 Effect of penalty factor on physical optical flow algorithm

實(shí)驗(yàn)3 圖4依次是基于物理光流算法和基于物理的加入各向異性擴(kuò)散后的光流算法得到的光流圖像，其中，擴(kuò)散系數(shù)方程參數(shù)α=5,β=0.8。由圖4可得，后者可顯著降低光流在邊緣地區(qū)的擴(kuò)散，光流較為稠密，但由于噪聲和異常點(diǎn)影響，光流穩(wěn)健性較差。

圖4 各向異性擴(kuò)散對(duì)基于物理光流算法的影響Fig.4 Effect of anisotropic diffusion on physical optical flow algorithm

實(shí)驗(yàn)4 圖5是使用基于物理光流算法和基于物理的穩(wěn)健光流算法得到的光流圖像，其中，α=5,β1=0.8,β2=0.02。由圖5可得，基于物理的穩(wěn)健光流算法可明顯減少邊緣光流擴(kuò)散，且顯著降低噪聲和異常點(diǎn)對(duì)光流計(jì)算的影響，得到較稠密光流，從而可改善光流計(jì)算的穩(wěn)健性。

圖5 基于物理光流算法和基于物理的穩(wěn)健光流算法對(duì)比

綜上所述，與LK,HS,基于金字塔LK以及基于物理光流算法相比，本文所提的引入各向異性濾波器及懲罰因子的基于物理的穩(wěn)健光流算法可增強(qiáng)光流計(jì)算對(duì)噪聲和異常點(diǎn)的穩(wěn)健性，同時(shí)明顯減少邊緣地區(qū)光流擴(kuò)散，進(jìn)而得到稠密的光流，從而可顯著改善光流計(jì)算的穩(wěn)健性。

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

針對(duì)可視化流動(dòng)圖像邊緣易擴(kuò)散和噪聲點(diǎn)及異常點(diǎn)的影響使得光流計(jì)算穩(wěn)健性較差的問(wèn)題，本文提出一種基于物理學(xué)的穩(wěn)健光流算法以改善光流計(jì)算的穩(wěn)健性。所提算法在基于物理光流算法中引入各向異性濾波器以增強(qiáng)邊緣光流的穩(wěn)健性，增加懲罰因子以減少噪聲及異常點(diǎn)對(duì)光流計(jì)算的影響。而后基于變分方法極小化光流能量函數(shù)求解歐拉-拉格朗日方程，最后通過(guò)迭代方法求得速度場(chǎng)。仿真結(jié)果表明，與傳統(tǒng)的LK,HS,基于金字塔LK以及基于物理光流算法相比，所提算法可顯著減少邊緣和角落區(qū)域的光流擴(kuò)散，改善針對(duì)噪聲及異常點(diǎn)的穩(wěn)健性，從而得到具有較好穩(wěn)健性的速度場(chǎng)。