基于光流算法的粒子圖像測速技術研究

王宏偉, 黃 湛

(中國航天空氣動力技術研究院, 北京 100074)

基于光流算法的粒子圖像測速技術研究

王宏偉, 黃 湛*

(中國航天空氣動力技術研究院, 北京 100074)

光流測量技術作為一種新的空氣動力學實驗技術，以其像素級分辨率的矢量場測量優勢獲得廣泛的應用。光流測量技術使用光流約束方程，配合平滑限定條件，可以進行速度場測量，獲得高分辨率的全局矢量場。首先通過研究積分最小化光流測速理論和算法，采用C++編寫光流速度測量程序，然后通過3種典型人工位移圖像對光流計算程序進行驗證，并將結果和標準位移分布進行比對分析，以指導如何在實際應用中獲得高精度光流速度場，最后進行小型風洞后向臺階實驗，利用高速相機拍攝示蹤粒子圖像，使用光流計算程序獲得速度矢量場，同采用互相關算法的粒子圖像測速計算結果進行比較，體現出光流計算方法像素級分辨率的矢量場測量優勢。

光流；速度測量；粒子圖像；運動估計；積分最小化

0 引 言

在空氣動力學發展過程中，應用實驗手段獲得流體的空間運動結構，獲取全場流動信息，是了解空氣動力學諸多流動現象及機理的關鍵，同時這一需求也不斷刺激空氣動力學實驗測試技術水平的提升。傳統的流動測量方法大多采用單點測量或介入式測量手段，不能滿足復雜流動情況以及全場速度測量的需求。

近年來，粒子圖像測速技術利用拍攝跟隨流場運動的粒子圖像并通過相關計算獲得速度場，并以其非介入式、全場流動速度測量、測量精度高及測速范圍廣等優點在全局速度測量方面具有應用優勢。不過，粒子圖像測速得到的速度矢量為其判讀區域內粒群的平均速度，忽略了此判讀區域內速度場的變化，這使得粒子圖像測速方法的空間分辨率受到一定的限制，尤其對于速度梯度大、速度變化劇烈的區域，其速度場測量精度是受到制約的。

相機所記錄的物體的圖像，本質上是物體所反射的光的分布。當相機連續拍攝運動物體時，物體所反射的光就會在相機底片上形成一系列連續變化的圖像，由這些圖像所組成的連續變化的信息不斷“流過”相機底片，如同光在流動，這是光流定義的物理基礎。物體在光的照射下，其表面灰度會呈現出一定的分布，而所謂光流，就是指圖像中灰度模式的運動速度。粒子圖像測速技術中拍攝到的粒子圖像顯示的是粒子濃度的分布，具有一定濃度的運動粒子群通過其散射光在相機底片上成像，所以粒子圖像測速技術中采用相關方法得到的速度場實際上是光流速度場(Optical Flow)，光流速度場能否代表真實速度場要取決于示蹤物質能否充分地跟隨流體運動。

Horn和Schunck早在20世紀80年代就開創性地進行了光流研究和計算[3]，Horn等認為同一運動物體引起的光流場應該是連續的、平滑的，即同一物體上相鄰點的速度是相似的，那么其投影到圖像上的光流變化也應該是平滑的，并提出了一種利用加在光流場上的附加速度平滑約束，即全局平滑約束將光流場的計算問題轉化為一個變分問題，給出像素級分辨率的運動估計。隨后，光流技術在發展過程中吸引了各個領域的科研人員進行研究[4-6]，也產生了很多其它解算方法，如圖像相關測速法、基于能量的方法[7]、基于相位的方法[8]和神經動力學方法等，針對Horn等提出的微分法在計算較大像素位移時遇到的困難，Tianshu Liu等在其工作中提到[9]，Corpetti等人提出在時間步上采用積分形式的光流方程代替微分形式的變分方程[10]，Ruhnau等使用高斯濾鏡生成二元圖像金字塔的多分辨率策略[11]。這些算法有各自的優勢，不過由于具有一定的適用范圍和算法復雜度使其相對于微分法無法獲得廣泛應用。

本文引入光流計算方法，編寫光流速度測量程序，通過人工位移圖像對該程序進行驗證和分析，以合理的測速范圍和最優的加權值指導如何在實際應用中獲得高精度光流速度場；進行小型風洞實驗，以粒子圖像為速度測量載體，使用光流計算方法獲得像素級分辨率的速度矢量場。旨在通過這些研究，發展以光流技術為基礎的全場流動測量技術，期望在提高空間分辨率和避免速度梯度影響、了解復雜流動的空間結構、獲取更多的動力學信息方面獲得進步，以此豐富實驗流體力學測量技術。

1 基于粒子圖像的光流速度場求解

Horn的光流場解算方法是微分法的基礎，它以灰度守恒方程為基礎，附加速度平滑約束條件，利用連續變化的圖像上各點灰度的時間、空間梯度來計算每一點的速度矢量。

光流的灰度約束條件認為，圖像上各點灰度隨時間和空間變化滿足如下關系：

式中：I(x,y,t)是t時刻(x,y)位置的圖像灰度，I(x+Δx,y+Δy,t+Δt)為原(x,y)位置的物體運動到(x+Δx,y+Δy,t+Δt)位置的灰度，令u=Δx/Δt,v=Δy/Δt，將上式進行泰勒展開并忽略高階項可以獲得光流約束方程：

式中：Ix,Iy,It是灰度的二維空間和時間梯度。顯然該方程是求解u,v2個未知數的不定方程，需要引入輔助約束方程進行求解。

速度梯度平方和判定量為

實際拍攝的圖像灰度將會受到硬件條件和噪聲的影響而偏離基本的光流約束條件，所以不能肯定地認為ξI恒為0，因此有必要找到合適的加權因子β，使得2個判定誤差在二維空間上的積分為最小，即對

ξ2=?

進行最小化。通過變分法建立每個因變量的歐拉特征方程：

2 光流速度場求解算法驗證

基于微分法編寫了光流速度場計算程序，擬用于風洞試驗速度場測量。為了驗證光流算法的精度，本文選取灰度圖像，通過人工位移獲得標準速度場，然后采用所編寫的光流計算程序對該圖像對進行處理，將光流計算結果與標準速度場進行比較。實際計算中，將采用1個單位的時間梯度，所得的速度值將直接反映圖像像素位移值，同時，也可考察2張圖像的像素位移對光流計算精度的影響。

首先考察第1組圖像(見圖1)，圖像為600pixel×600pixel大小的風暴云圖，通過處理將圖像向右平移一個像素(Δx=1pixel)，通過光流速度場計算程序得到的速度矢量分布(Skip 16pixel×16pixel，下同)和流線如圖2所示，可以看出，在1像素位移條件下，光流計算方法所得到的速度場分布十分均勻。

圖1 灰度原圖與平移一個像素后圖像

圖2 采用光流計算方法得到的速度矢量圖和流線圖

圖3 橫向速度分布和縱向速度分布分析

在圖1中所示的縱向考察線和橫向考察線(下同)上對光流速度場和標準速度進行分析如圖3所示，可以看出, 在1像素位移條件下光流算法獲得的速度值與標準速度幾乎沒有偏差。

第2組圖像(見圖4)采取了Δx=1 pixel和Δy=1 pixel的斜向移動，處理得到的速度矢量分布和流線分布如圖5所示，速度矢量分布很均勻，流場平滑性很好，不過已開始出現小范圍的速度偏差。

通過縱向考察線和橫向考察線將所計算的速度場與標準速度進行比較發現(見圖6), 在Δx=1 pixel和Δy=1 pixel的斜向運動條件下，所計算的速度值能夠保證較高的準度分布，不過已開始出現一定量的抖動偏差。這種偏差很有可能與2個方向速度偏差的合成有關；同時也暗示隨著連續圖像間運動位移的增大，光流算法的計算誤差也會隨之增加。

圖4 灰度原圖與x,y各平移一個像素后圖像

圖5 采用光流計算方法得到的速度矢量圖和流線圖

圖6 橫向速度分布和縱向速度分布與標準速度比較

第3組圖像(見圖7)為中心旋轉1°前后的2張圖像，通過計算獲得的速度矢量和流線分布圖如圖8所示，可以看出，圖像中心的速度矢量分布較為均勻和平滑，邊緣處的速度偏差較大，平滑性較差。

通過在前述2個方向上將光流計算的速度矢量和標準速度進行比較(見圖9)，可以明顯發現：圖像位移量不超過2個像素部分的計算結果與標準速度結果十分吻合，偏差很?。粓D像位移量超過2個像素部分的計算結果則出現了明顯的偏差。這很好地說明了所得到的速度場分布中為何中心區域速度分布均勻，平滑性好，而邊緣區域速度偏差較大，平滑性較差，同時也在很大程度上對如何在實際測量應用中保證光流計算獲得高精度的結果進行了指導。

圖7 灰度原圖與旋轉1度后的圖像

圖8 采用光流計算方法得到的速度矢量圖和流線圖

圖9 橫向速度分布和縱向速度分布與標準速度比較

3 加權因子最優值估計

β2作為加權因子表征了實驗和數值離散噪音的相對大小。在實際計算中β2通過調整速度梯度平方和比重來影響計算結果的平滑性，即通過選取不同的加權值獲得不同平滑程度的流場速度分布。下面將以旋轉灰度圖像為計算基礎，在不同β2值下對光流計算結果中縱向考察線速度分布與標準速度分布進行比較(見圖10～14)。

由圖中數據可以看出，隨β2增大，縱向考察線上速度場的分布能夠更好地吻合標準速度分布。在實際應用粒子圖像進行光流速度場計算時，β2相對較大，說明噪聲較大，這與粒子圖像不同于激光誘導熒光這類圖像有關，粒子圖像的灰度分布存在一定的間斷性，需要大的加權值來限定誤差。

圖10 β2=1時縱向考察線速度分布與標準速度分布比較

圖11 β2=4時縱向考察線速度分布與標準速度分布比較

圖12 β2=7時縱向考察線速度分布與標準速度分布比較

圖13 β2=10時縱向考察線速度分布與標準速度分布比較

圖14 β2=13時縱向考察線速度分布與標準速度分布比較

為了進一步分析β2的選取對速度場特性影響，根據不同β2值解算出的光流速度場與標準速度分布之間的相關系數來表現β2對解吻合的影響(見圖15)。根據圖中數據結果可以看出，解的準確性在β2增長初期會迅速增加，在10～18之間達到平穩階段并取得極值，超過20以后，反而會下降，因此選取適當的β2值將對光流算法獲得的速度場的準確性有很大幫助。不過這里β2是與灰度的時間和空間梯度間隔相關的常數，實際應用中應結合具體參數進行選取計算。

圖15 相關系數隨β2變化曲線

4 光流速度場測量試驗與結果

本文針對光流測速算法對流場的要求，進行了后向臺階內流場渦結構觀測實驗。

4.1 光流速度場測量試驗

實驗的難點在于如何保證實驗條件滿足光流流場計算的約束方程。根據前述分析，應用如上微分法對光流流場進行解算時，要求連續拍攝的每張圖像具有相同的光照條件，相鄰圖像對應像素點的位移不至太大，最好保證在2個像素以內，從而進行有效的光流速度場計算。

實驗在實驗室自行搭建的小型低速風洞系統內進行，整個系統由風洞主體、流動控制系統、示蹤粒子發生和播發系統、激光照明及附屬光路系統、圖像采集系統和計算機處理系統組成。實驗布局如圖16所示。

圖16 光流測速實驗布局

風洞主體通過其前段可控轉速的渦輪風扇電機控制流動速度，通過調節變頻器可以將流經中部實驗段的流動速度控制在1～20m/s的范圍內，風洞主體通過其尾部的蜂窩整流器吸入混有示蹤粒子的氣流并進行整流，保證來流的穩定條件，并與3個示蹤粒子播撒管共同工作使示蹤粒子與吸入氣流混合均勻。為保證恒定的光照條件，實驗采用氬離子激光器輸出的連續激光通過擴束鏡和片光源為流場提供連續的片光照明；為保證相鄰圖像間位移不至太大，實驗采用高速相機進行圖像采集，觀測流場前先對標尺進行拍攝用以標定空間尺度，采集結束后傳輸到計算機上存儲并處理,現場設備及布置照片如圖17～20所示。

圖17 實驗光路及實驗段

圖18 粒子播撒架

圖19 粒子發生器

4.2 粒子圖像光流場分析與圖像均一化

高速相機拍攝的后向臺階內流場粒子圖像如圖21所示，大量連續拍攝的粒子圖像記錄了連續變化的光流流場信息，虛線框內為計算區域。為了對實際拍攝示蹤粒子圖像的光流場特性進行驗證，將實際拍攝的光流場圖像進行旋轉獲得標準位移場，如圖22所示；采用與前述光流解算方法驗證相同的比較方式對計算獲得的速度場進行考察，如圖23和24所示；并分析了速度分布與標準速度分布間的相對誤差和絕對誤差，如圖25所示。

圖20 實驗現場

圖21 后向臺階流場粒子圖像及流場計算區域

圖22 使用示蹤粒子圖像旋轉1°前后圖像及驗證計算區域

圖23 用光流算法對旋轉前后粒子圖像進行計算的速度矢量與流線圖

分析可以看出，在示蹤粒子圖像上，當像素位移小于一定值時，計算獲得的速度分布能夠很好的與標準值吻合；橫向考察線和縱向考察線上速度分布絕對誤差基本處于5%線以下，不過由于靠近0位移位置時，速度絕對值較小，故而相對誤差較大；在±1像素位置絕對誤差和相對誤差最小，都達到了1%左右，該位置±0.5像素鄰域是速度相對精度和絕對精度都比較高的區域。

圖24 橫向速度分布和縱向速度分布與標準速度比較

圖25 考察線上速度分布的絕對誤差和相對誤差

由于激光經過柱透鏡制造的片光存在一定的擴張角度，在激光片光傳播過程中，片光的寬度始終在變大，所以激光片光在展開區域的能量分布隨時空變化也會有一定的差異。此外由于實驗段的有機玻璃制造工藝以及存在一定的磨損，其中存在許多細微的氣泡或表面有輕微劃痕，當激光入射進有機玻璃遇到這些氣泡或劃痕時，激光就要被削減一部分能量，導致射入流場里的激光片光含有大量的條紋，即激光片光的能量分布是不連續的，自然所得的灰度場里也含有大量的條紋。

如果還要考慮噪聲的話，可以在實驗前或實驗后獲得一個噪聲灰度場。獲得噪聲灰度場的一種方法是，在不加示蹤物質的情況下，將攝像機對準拍攝區域，連續拍攝一段時間，將得到的圖像進行平均，就可以得到噪聲灰度場INoise。這時相對灰度場的表達如式(9)所示。

4.3 試驗結果與分析

經過光流算法計算得到的真實流場速度分布如圖26所示，結果顯示最大位移為2.2像素，表明區域內速度計算結果基本位于可信區間內。

由于所觀測流場的絕對速度較低(<1m/s),后向臺階內流場結構不具有典型性，所以將計算區域內的光流處理結果同采用互相關算法的PIV結果(見圖27)進行比較。圖28為PIV顯示的流場結構，同圖29的光流算法流場結構比較可以發現，在同樣像素分辨率下，光流算法的結果在速度大小、方向和等速度線分布上基本與PIV計算結果基本一致，大尺度渦結構和位置也與其基本無差異，但光流計算結果中可以觀測到更小尺度的渦結構，等速度線更加平滑。

圖26 光流計算結果(Skip 8pixel×8pixel)

圖27 PIV計算結果(8pixel×8pixel判讀區)

比較兩者流線圖(圖28和29)可以看出，在同樣像素分辨率下，PIV算法的計算結果中較大尺度的渦結構可以顯示出來，但是對于比較小尺度的渦結構顯示卻不明顯，而光流算法的計算結果中小尺度渦結構被明顯地顯示出來。

為進一步分析，將兩者小尺度渦結構區域圖像放大(圖30和31)，可以看出，光流計算結果在每個像素上都獲得一個速度矢量，可以顯示更加精細的流場結構，而PIV算法由于采用了8pixel×8pixel的判讀小區，其流場結構的精細顯示受到限制。

圖28 PIV計算獲得的流場結構

圖29 光流計算獲得的流場結構

圖30 PIV結果中小尺度結構區域放大

圖31 光流結果中小尺度結構區域放大

5 結論與展望

通過對光流算法研究與編程實現，成功搭建了可用于光流測速的實驗平臺，進行了小型風洞內的光流速度場測量實驗，比較光流算法和PIV對拍攝的粒子圖像處理得到的速度場，可以得出結論：在本實驗流動條件下，光流算法可以在像素級分辨率情況下獲得優于PIV獲得的速度矢量場，適合用于復雜流動的速度場測量。不過基于微分法的光流解算方法在計算較大位移區間時存在缺陷，未來將進一步通過算法研究和實驗手段進行改進。

[1] Lang D B. Laser doppler velocity and vorticity measure- ments in turbulent shear layers[D]. California:Caltech, Pasadena, 1985.

[2] Adrian R J. Particle-imaging techniques for experimental mechanics[J]. Annual Review of Fluid Mechanics, 1991, 23(1): 261-304.

[3] Horn B K P, Schunck B G. Determining optical flow[J]. Artificial Intelligence, 1981, 17: 185-203.

[4] Liu T, Shen L. Determination of velocity and skin friction fields from images by solving projected motion equations[C]. 22nd International Congress on Instrumenta- tion in Aerospace Simulation Facilities (ICIASF), International Congress on Instrumentation in Aerospace Simulation Facilities, Pacific Grove, CA, June 2007.

[5] Liu T, Sullivan J P. Pressure and temperature sensitive paints experimental fluid mechanics[M]. Berlin: Springer-Verlag, 2004.

[6] Tikhonov A N, Arsenin V Y. Solutions of ill-posed problems[M]. Winston, Washington D C, 1977.

[7] Heeger D J. Model for the extraction of image flow[J]. Optical Society of America, 1987, A4: 1455-1471.

[8] Fleet D J, Jespon A D. Computation of component image velocity from local phase information[J]. International Journal of Computer Vision, 1990, 5: 77- 104.

[9] Liu Tianshu, Shen Lixin. Fluid flow and optical flow[J]. J. Fluid Mech, 2008, 614: 253-291.

[10] Corpetti T, Heitz D, Arroyo G, et al. Fluid experimental flow estimation based on an optical flow scheme[J]. Exps Fluids, 2006, 40(1): 80-97.

[11] Ruhnau P, Kohlberger T, Schnorr C, et al. Variational optical flow estimationfor particle image velocimetry[J]. Exps Fluids, 2005, 38(1): 21-32.

[12] Zhan Huang, Hong Wei Wang, et al. The research of particle image velocimetry based on optical flow[C]. 16th International Symposium on Flow Visualization, Okinawa, Japan, 2004.

(編輯：楊 娟)

Research on particle image velocimetry based on optical flow

Wang Hongwei, Huang Zhan*
(China Academy of Aerospace Aerodynamics, Beijing 100074, China)

As a new aerodynamics experimental technique, optical flow test technique gradually attracts more and more attention due to its advantages in vector field measurement with pixel scale resolution and strong smoothness ability. By means of scalar constraint equation combined with smoothness constraint condition, optical flow test technique can measure global velocity vector field with high space resolution. In this paper, the integration minimization optical flow velocimetry theory and algorithm were firstly studied and programmed by C++; then, in order to verify the accuracy of the optical flow algorithm, three groups of grayscale images shifted by given displacements were used for processing by the optical flow algorithm program and the result can be used to guide how to obtain high-precision optical flow calculation result in the actual measurement applications; at last, backward step verification experiment, in which tracer particle images were acquired by high speed camera and velocity vector field was calculated by optical flow algorithm, was completed in a small wind tunnel. The result calculated by optical flow algorithm was compared with the result calculated by PIV′s cross-correlation algorithm, which shows that optical flow test technique possesses the advantages of pixel scale resolution.

optical flow;velcocity measurement;particle image;motion estimation;integration minimization

1672-9897(2015)03-0068-08

10.11729/syltlx20140115

2014-10-11;

2014-12-02

國家重點基礎研究發展計劃(2014CB744801)

WangHW,HuangZ.Researchonparticleimagevelocimetrybasedonopticalflow.JournalofExperimentsinFluidMechanics, 2015, 29(3): 68-75. 王宏偉, 黃 湛. 基于光流算法的粒子圖像測速技術研究. 實驗流體力學, 2015, 29(3): 68-75.

O355

A

王宏偉(1987-)，男，黑龍江訥河人，助理工程師。研究方向：實驗流體力學。通信地址：北京市豐臺區云崗西路17號(100074)。E-mail: www_whw_cn@163.com

*通信作者 E-mail: xunfang05@sohu.com