鄭 佳， 王洪雁， 裴炳南

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

0 引言

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

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

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

1 經(jīng)典光流法

1.1 光流約束方程

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

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

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

(1)

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

(2)

Ixu+Iyv+It=0

(3)

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

1.2 Lucas-Kanade光流法

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

(4)

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

A2W2Ap=ATW2b

(5)

式中，

(6)

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

p=(ATW2A)-1ATW2b。

(7)

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

1.3 Horn-Schunck光流法

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

(8)

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

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

(9)

基于上述兩個約束條件可知，光流矢量可通過如下最小化問題獲得，即

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

(10)

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

2 基于物理的光流法

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

(11)

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

(12)

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

(13)

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

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

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

(14)

(15)

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

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

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

(16)

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

(17)

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

加入懲罰項，式(17)可改寫為

(18)

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

(19)

式中：

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

(20)

(21)

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

(22)

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

(23)

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

(24)

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

對式(24)整理可得

(25)

式中，

(26)

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

(27)

式中，

(28)

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

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

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

3) 采用九點差分格式計算u，v均值[15]。

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

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

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

4 實驗仿真及分析

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

(29)

式中，

(30)

實驗1 圖1所示為等流速狀態(tài)下兩幀連續(xù)的粒子圖像，以下算法皆基于此計算光流。圖2所示為LK光流算法、HS光流算法、基于金字塔LK光流算法以及基于物理算法得到的光流圖像。由圖2可知，與HS光流算法相比，LK光流算法得到的光流稀疏，穩(wěn)健性較好。HS光流算法雖可得稠密光流，但穩(wěn)健性差。較HS光流算法和LK光流算法，基于金字塔LK光流算法所得光流精度較高，穩(wěn)健性較強，但光流較稀疏。此3種光流算法所得光流分辨率較低。基于物理光流算法可得高分辨光流，但邊緣易擴散，且易受噪聲和異常點影響，光流較稀疏，穩(wěn)健性較差。

圖1 原始圖像Fig.1 Original images

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

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

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

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

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

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

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

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

5 結(jié)束語

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