李 超蔣暑民馬洪余戴德君黃傳江

(1.自然資源部 第一海洋研究所,山東 青島266061;

2.青島海洋科學(xué)與技術(shù)試點國家實驗室 區(qū)域海洋動力學(xué)與數(shù)值模擬功能實驗室,山東 青島266237;3.自然資源部 海洋環(huán)境科學(xué)與數(shù)值模擬重點實驗室,山東 青島266061)

激光粒子圖像測速技術(shù)(Particle Image Velocimetry,PIV)是采用高速CCD相機、激光器等設(shè)備對流場進行成像,再通過對相鄰兩幅流場圖像加以互相關(guān)分析獲取測量區(qū)域的流速和流向的方法[1]。相較于傳統(tǒng)的單點測量方法,PIV技術(shù)具有多點同時測量、非接觸測量、瞬時性和高分辨率等優(yōu)點。因此,自20世紀80年代PIV技術(shù)被Adrian等學(xué)者提出以后迅速發(fā)展[2-5],并被應(yīng)用在多個方面,如燃燒的火焰場、內(nèi)燃機、微血管、飛機機翼復(fù)雜的外形流動以及火炮發(fā)射口流場測量等。

在應(yīng)用PIV過程中,首先,需要向待測流場中均勻散播示蹤粒子,這些粒子應(yīng)具備如下特性:反光性良好、跟隨性良好、不改變流體性質(zhì),各個示蹤粒子之間互不影響等[6-8];然后,利用脈沖激光或連續(xù)激光片光源照亮待測量的流場區(qū)域,高速CCD相機攝取激光面上的粒子分布,在極短的時間內(nèi)記錄粒子圖像幀序列;最后,選定合適的窗口針對粒子圖像序列進行互相關(guān)分析和亞像素算法分析,計算得出2幀圖像中相對應(yīng)的粒子位移的大小和方向,從而獲得測量區(qū)域的瞬時速度場[9-10]。在這個過程中,圖像匹配是整個測量過程中的關(guān)鍵技術(shù)之一,示蹤粒子直徑的大小和粒子的濃度對整個測量結(jié)果也起著至關(guān)重要的作用[11]。

相鄰圖像的相關(guān)分析或圖像匹配是PIV技術(shù)獲取流場信息的關(guān)鍵,對原始圖像的相關(guān)分析[12-13]通常分為2步:1)通過相關(guān)搜索匹配獲得整像素位移,2)在整像素位移基礎(chǔ)上進行亞像素的位移匹配。整像素位移很容易獲得,但在實際流場中,所選匹配區(qū)域的位移值一般不會恰好為像素的整數(shù)倍,所以整像素的位移精度在實際應(yīng)用中遠遠不夠[14]。為了提高PIV的測量精度,通常采取3種方法:1)提高高速CCD相機的分辨率,但是由于當前相機技術(shù)限制,CCD相機的分辨率有限,很多時候并不能滿足需求;2)提高攝像系統(tǒng)的放大倍數(shù),但放大倍數(shù)的增加意味著可測量面積的減小;3)采用亞像素位移匹配算法,亞像素位移算法在高精度PIV測量中具有重要的意義[15-16]。

本文比較了二次多項式曲面擬合法和基于梯度的亞像素位移匹配算法,給出2種算法的匹配精度和各自的特點,并進一步評估了匹配窗口的大小、示蹤粒子的直徑和粒子的濃度對測量精度的影響。相關(guān)結(jié)果對于合理使用PIV技術(shù)開展實驗和測量結(jié)果誤差分析具有一定的參考價值。

1 數(shù)據(jù)與方法

1.1 計算機仿真粒子圖像

為了比較不同的亞像素定位算法精度且找到合適的粒子大小和粒子濃度,必須采用能夠精確控制位移的粒子圖像。本文采用Zhou和Goodson[17]提出的算法,利用計算機產(chǎn)生隨機分布的高斯光斑來模擬粒子的分布,并精確控制粒子的移動距離。仿真粒子圖像的生成函數(shù)可表示為

式中:x,y為仿真粒子圖像每個像素灰度值所對應(yīng)的位置;I1,I2分別是原始仿真粒子圖像和精確移動后的仿真粒子圖像灰度值;△u,△v為原始圖像仿真粒子移動的位移分量;s為仿真粒子的數(shù)目;I0為仿真粒子的中心光強;x k,y k為仿真粒子在原始圖像中的中心位置;a為仿真粒子的直徑。

利用式(1)、式(2)生成2幅仿真粒子圖像,如圖1所示:其中,仿真粒子圖像大小為128像素×128像素;仿真粒子數(shù)目為300個,中心光強為255,直徑為4個像素。圖1a為原始仿真粒子圖像,圖1b是在圖1a基礎(chǔ)上所有粒子向右移動1.5個像素后的仿真圖像,即△u=1.5且△v=0的情況。

圖1 仿真粒子圖像Fig.1 Particle images simulated by computer

1.2 亞像素位移匹配算法

1.2.1 互相關(guān)匹配算法獲取整像素位移

亞像素位移的獲取首先通過互相關(guān)匹配獲得整像素位移,然后再基于整像素匹配結(jié)果進行亞像素匹配?；ハ嚓P(guān)匹配需要2幅粒子圖像,在前一時刻粒子圖像中選取匹配窗口f(x,y),在后一時刻粒子圖像中選取同樣大小的匹配窗口g(i,j)進行匹配,其中相關(guān)系數(shù)最大的窗口位置即可視為2幀照片的整像素位移量。本文整像素位移的獲取采用Willert和Gharib[18]在1991年提出的相關(guān)函數(shù):

式中,R(i,j)為相關(guān)系數(shù);(x,y)為前一時刻匹配窗口內(nèi)灰度值對應(yīng)位置坐標;i,j為后一時刻匹配窗口相對于前一時刻窗口的整像素移動位置,一般選擇窗口大小為32像素×32像素或64像素×64像素;通過上述相關(guān)函數(shù)獲取整像素位移之后,需要進一步通過亞像素位移匹配算法來提高測量精度。本文比較了曲面擬合法和基于梯度的亞像素位移算法。

1.2.2 曲面擬合法

實際計算中常用的曲面擬合法為二次多項式曲面擬合,擬合的參量為獲取整像素位移時計算得到的相關(guān)系數(shù)。二次多項式曲面擬合法需要先選定擬合窗口,通常選取擬合窗口為3×3,即根據(jù)整像素位移結(jié)果及其周圍相近8個點的相關(guān)系數(shù)組成3×3的相關(guān)系數(shù)矩陣,利用二次多項式擬合為連續(xù)曲面,找出此曲面的極值點所對應(yīng)的位置即為最佳匹配結(jié)果。

二次多項式的解析方程為

式中,C為曲面擬合相關(guān)系數(shù),(x,y)為以整像素位移匹配結(jié)果為中心點建立坐標系中的坐標點,用最小二乘法求解待定系數(shù)k i,進而可求得極大值點,該極大值點對應(yīng)位置即曲面擬合法獲得的亞像素精度的位移值。

1.2.3 基于梯度的亞像素位移算法

鑒于在PIV測量過程中,粒子位移量不大,假設(shè)位移前后同一點的灰度值相同,則:

式中,f(x,y)和g(x,y)分別為前一時刻和后一時刻所對應(yīng)的匹配窗口內(nèi)的灰度值;△x,△y為整像素位移值;x0和y0為亞像素位移值。將式(5)進行一階泰勒展開并舍去高階小量,可得:

假設(shè)匹配窗口m×m內(nèi)各個像素點的位移都相等,就有m×m個式(6),再利用最小二乘法可以推導(dǎo)出:

求得亞像素位移值x0和y0并結(jié)合整像素位移即可獲得精度較高的位移值。

2 結(jié)果與分析

2.1 曲面擬合法與基于梯度亞像素位移算法精度對比

前文介紹了曲面擬合法和基于梯度亞像素位移算法,下面將對這2種亞像素算法的匹配精度和特征進行對比。首先針對粒子均勻移動的理想情況進行分析,所謂均勻移動指的是圖像中的所有粒子移動大小和方向完全一致。利用的仿真粒子圖像大小為128像素×128像素,仿真粒子數(shù)目為1000個,仿真粒子直徑為3個像素。圖像互相關(guān)匹配時的匹配窗口大小為32像素×32像素。圖像互相關(guān)匹配前進行簡單的濾波,并不會影響仿真粒子的直徑大小。圖2給出了二次多項式曲面擬合法與基于梯度亞像素位移算法的匹配精度及特征。

圖2 二次多項式曲面擬合法與基于梯度亞像素位移算法得到的亞像素匹配結(jié)果及誤差對比Fig.2 Subpixel displacements and errors by using quadratic polynomial surface fitting and subpixel displacement algorithm based on gray gradient

圖2a和圖2b為使用二次多項式曲面擬合法得到的結(jié)果,圖2c和圖2d為使用基于梯度的亞像素位移算法得到的結(jié)果。針對選定的真實位移量,通過選取不同的圖像對和在前一時刻圖像中選取不同窗口位置,在下一時刻圖像中進行匹配,獲得多個匹配結(jié)果,進而進行誤差分析。圖2為針對每個選定的真實位移量,100次匹配實驗結(jié)果:二次多項式曲面擬合法整像素位移時的誤差較小,約為0.02個像素,而在非整像素位移時誤差比較大(圖2a),且平均誤差值隨著粒子真實位移的改變呈現(xiàn)出周期性的變化(圖2b),但均在0.08個像素內(nèi);基于梯度亞像素位移算法在整像素時的誤差為0,所有誤差也在0.08個像素左右(圖2c,圖2d)。

多次匹配實驗結(jié)果的平均絕對誤差分析(表1)可以看出:二次多項式曲面擬合法和基于梯度亞像素位移算法在相同匹配窗口時x方向和y方向的匹配精度差別不大。二次多項式曲面擬合法的計算精度隨匹配窗口的增大而提高,基于梯度亞像素位移算法在不同匹配窗口下的平均絕對誤差大致相同。在匹配窗口為16像素×16像素時,基于梯度亞像素位移算法的匹配精度要優(yōu)于二次多項式曲面擬合法的匹配精度;匹配窗口為32像素×32像素時,2種算法的匹配精度大致相同;匹配窗口為64像素×64像素時,二次多項式曲面擬合法要優(yōu)于基于梯度亞像素位移算法。

表1 不同窗口大小下兩種亞像素匹配平均絕對誤差比較Table 1 Mean absolute errors of subpixel displacements by using quadratic polynomial surface fitting and subpixel displacement algorithm based on gray gradient for different window sizes

2.2 粒子濃度對測量結(jié)果的影響

在PIV測量流體過程中需要向流體中撒播示蹤粒子,粒子濃度對PIV的測量精度有著十分重要的影響。利用準確控制仿真粒子圖像的位移和不斷改變生成粒子的數(shù)目,可以找到合適的粒子濃度范圍,對實際實驗過程中的PIV測量有一定的指導(dǎo)意義。圖3給出了采用二次多項式曲面擬合法時不同粒子濃度對匹配精度的影響,采用的仿真粒子圖像大小為128像素×128像素,仿真粒子直徑為3個像素,仿真粒子真實位移為向右和向上分別移動1.5個像素。選擇合適的匹配窗口大小,通過選取不同的圖像對和在前一時刻圖像中選取不同窗口位置,在下一時刻圖像中進行匹配,獲得100個匹配結(jié)果,進而統(tǒng)計得到誤差分析結(jié)果。

圖3給出了匹配窗口為16像素×16像素時,匹配誤差隨仿真粒子濃度的變化,這里粒子濃度指仿真粒子個數(shù)與仿真粒子圖像總像素值的比值。從圖中可以看出:匹配窗口為16像素×16像素時,粒子濃度0.04～0.36最為合適。改變匹配窗口大小后,最合適的仿真粒子濃度區(qū)間也發(fā)生改變,匹配窗口為32像素×32像素時,最合適的粒子濃度區(qū)間為0.008～0.440,匹配窗口為64像素×64像素時,最合適的粒子濃度區(qū)間為0.002～0.530。也就是說,匹配窗口越大,其包含的灰度值信息越多,最合適的粒子濃度范圍也越大,建議在實際實驗時,同時考慮所研究科學(xué)問題對流場空間分辨率的要求(匹配窗口的大小)和實驗經(jīng)濟性,選擇合適的粒子濃度。

圖3 亞像素匹配誤差隨仿真粒子濃度變化圖(16像素×16像素)Fig.3 Variation of subpixel displacement errors with particle concentration(16 pixels×16 pixels)

需要說明的是,這里所指的粒子濃度是仿真粒子個數(shù)與仿真粒子圖像總像素值的比值。在利用計算機生成仿真圖像時,設(shè)置了避免粒子重合的算法,但也只是保證仿真粒子中心位置不會重合,導(dǎo)致在生成的仿真圖像中,會不可避免的產(chǎn)生粒子邊緣重合現(xiàn)象。因而本文給出的結(jié)果對于實驗中所需的最低粒子濃度具有一定的指導(dǎo)意義。

2.3 粒子直徑對測量結(jié)果的影響

PIV是一種非接觸式測量技術(shù),粒子的選擇對測量結(jié)果有比較大的影響,首先粒子的密度必須與所測量流體的密度接近,這樣選取的示蹤粒子才能具有較好的跟隨性,能很好地反映流場信息,粒子直徑大小對測量精度也有很大的影響。我們利用精確控制仿真粒子的大小并采用二次多項式曲面擬合法來分析測量誤差的變化,進而確定適合的粒子大小。仿真粒子圖像大小為128像素×128像素,仿真粒子的數(shù)目為1000個,對應(yīng)的粒子濃度為0.061,仿真粒子位移為向右和向上分別移動1.5個像素。通過選取不同的圖像對和在前一時刻圖像中選取不同窗口位置,在下一時刻圖像中進行匹配,獲得100個匹配結(jié)果,進而統(tǒng)計得到誤差。

圖4給出了不同匹配窗口下不同粒子直徑時的誤差大小。在匹配窗口為16像素×16像素時,粒子直徑為2個像素和3個像素時誤差較小,其中粒子直徑為3個像素時誤差最小。當匹配窗口為32像素×32像素時,粒子直徑合適的范圍為3～5個像素,在粒子直徑為3個像素時誤差最小,其誤差值在0.02個像素內(nèi)。當匹配窗口變?yōu)?4像素×64像素時,粒子直徑合適的范圍為3～6個像素,粒子直徑在5個像素時匹配精度最佳,其中粒子直徑為3,4,5個像素時的誤差值在0.02個像素內(nèi)?？傊?隨著匹配窗口的增大,粒子直徑合適的范圍也在擴大。從整體來看,在這3種匹配窗口下,粒子直徑為3個像素時的匹配誤差均可接受。一般來講,粒子直徑越小,其跟隨性越好,尤其是在測量湍流流場條件下,對粒子的跟隨性的要求更高。

圖4 不同匹配窗口下亞像素匹配誤差隨粒子直徑變化圖Fig.4 Variation of subpixel displacement errors with particle diameter for different window sizes

3 理想渦旋流場下的誤差分析

前文在亞像素匹配精度分析時均采用粒子均勻移動的圖像,匹配窗口內(nèi)部位移相同,其運動方式較為簡單,而在PIV的現(xiàn)實應(yīng)用中很少會出現(xiàn)這種理想流動情況。因此,采用存在理想渦旋的粒子圖像,分析其亞像素匹配誤差更具有實際應(yīng)用價值,下面對這種復(fù)雜的流場情況進行分析。仿真粒子圖像大小為256像素×256像素,仿真粒子數(shù)目為3500個,對應(yīng)的粒子濃度為0.053,仿真粒子直徑為3個像素,仿真粒子的位移采用線性增加的方法,即距離旋渦中心越遠,粒子的位移越大且方向垂直于旋渦中心指向該點的矢量。粒子圖像互相關(guān)匹配時的匹配窗口大小為16像素×16像素,圖5為二次多項式曲面擬合法與基于梯度亞像素位移算法得到的粒子影像及反演渦旋流場圖,可以看出二者的差別不大,都能很好的表征渦旋的流場情況。

圖6給出了使用二次多項式曲面擬合法和基于梯度亞像素位移算法得到誤差分析結(jié)果,由圖6可見:使用基于梯度亞像素位移算法得到的誤差比使用二次多項式曲面擬合法得到的誤差要小,使用二次多項式曲面擬合法的匹配誤差基本穩(wěn)定在0.2個像素以內(nèi),使用基于梯度亞像素位移算法得到的誤差基本在0.05個像素以內(nèi)。在粒子沿直線運動時,且匹配窗口為16像素×16像素情況下,基于梯度亞像素位移算法匹配精度高于二次多項式曲面擬合法。在復(fù)雜的流動狀態(tài)下,匹配窗口內(nèi)部存在粒子的相對位移時,基于梯度亞像素位移算法的匹配精度也明顯優(yōu)于二次多項式曲面擬合法。

圖5 二次多項式曲面擬合法與基于梯度亞像素位移算法得到的渦旋流場Fig.5 Particle image and the vortex flow field obtained by using quadratic polynomial surface fitting and subpixel displacement algorithm based on gray gradient

圖6 二次多項式曲面擬合法與基于梯度亞像素位移算法得到的誤差對比Fig.6 Errors by using quadratic polynomial surface fitting and subpixel displacement algorithm based on gray gradient for the idealized vortex flow

4 PIV在水槽實驗中的應(yīng)用

2012年自然資源部第一海洋研究所建設(shè)了風(fēng)-浪-流多功能實驗水槽,水槽長45.0 m,寬1.0 m,高1.8 m,實驗時水深一般設(shè)為1.2 m,為實現(xiàn)平穩(wěn)出流,水槽兩端各設(shè)增壓腔一個。利用該水槽開展了波湍相互作用等實驗研究。圖7為風(fēng)-浪-流多功能實驗水槽玻璃段的示意圖(玻璃段長32.4 m),造波機安裝在水槽的首端,消波板安裝在水槽尾部。在開展波-湍相互作用實驗時,利用格柵振動產(chǎn)生湍流,采用PIV和ADV(Acoustic Doppler Velocimetry)測量水體速度,激光片光源平行于水槽,振動格柵位于測量區(qū)域的下方[20]。

圖7 風(fēng)-浪-流多功能實驗水槽玻璃段示意圖Fig.7 Schematic diagram of wind-wave-flow multifunctional tank

圖8給出的是格柵在水面以下26～31 cm處上下振動(頻率為3 Hz,振幅為5 cm)時垂直采樣得到的PIV影像。PIV測量時,以激光脈沖作為光源,采用高速CCD相機記錄粒子圖像,每兩幀粒子圖像之間的時間間隔為5 ms,采樣頻率為14.5 Hz,即每兩對粒子圖像之間時間間隔為0.069 s,粒子圖像分辨率為1192像素×1600像素,對應(yīng)的物理空間約為9.5 cm×12.8 cm,PIV影像的上端距水面約為4 cm。以其中一對粒子圖像為例進行分析,對所選擇的兩幀粒子圖像進行簡單的濾波后識別出粒子直徑的平均大小為2.96個像素,粒子直徑不是整數(shù),主要原因在于所用的示蹤粒子在原始物理形態(tài)上就難以保證大小均勻,另外在利用圖像識別粒子大小時,采用先識別每個粒子的面積,再假設(shè)每個粒子為圓形計算得出的粒子直徑大小。PIV影像中粒子個數(shù)約為10600個,對應(yīng)粒子濃度為0.0056?；ハ嚓P(guān)匹配窗口選取64像素×64像素,利用二次多項式曲面擬合法和基于梯度亞像素位移算法進行亞像素匹配,結(jié)果如圖8所示。2種算法反演得到的速度場除了個別點有差別之外,大部分區(qū)域都幾乎一致。

圖8 水槽格柵振動產(chǎn)生湍流實驗PIV粒子影像及反演流場Fig.8 Particle image and the flow field of grid-generated turbulence in wave tank

我們選取了格柵振動試驗中測得的200對連續(xù)的粒子圖像,比較了2種算法所得速度場。圖9給出了某一點(距水面14.720 cm,距測量區(qū)域左側(cè)4.576 cm)2種亞像素算法所得速度的對比圖,u為水平方向流速,v為垂直方向流速?；ハ嚓P(guān)匹配窗口為64像素×64像素,兩種亞像素算法所得流速變化趨勢與量值基本一致,在使用合適的粒子大小和濃度時,2種亞像素算法均能得到較好的結(jié)果。

圖9 二次多項式曲面擬合法與基于梯度亞像素位移算法得到的速度對比Fig.9 Comparisons of velocities obtained by using quadratic polynomial surface fitting and subpixel displacement algorithm based on gray gradient

5 結(jié) 論

本研究利用計算機仿真PIV粒子圖像,通過精確控制粒子位移,討論了不同匹配窗口下二次多項式曲面擬合法與基于梯度亞像素位移算法的匹配精度;進一步討論了粒子濃度和直徑大小對PIV測量結(jié)果的影響,給出了合適的粒子濃度范圍和最佳的粒子直徑大小;進而分析了理想旋渦流場下2種亞像素算法的誤差,最后將2種算法應(yīng)用到實際水槽實驗PIV影像分析中。具體結(jié)論如下:

1)比較了二次多項式曲面擬合法與基于梯度亞像素位移算法2種方法,給出了不同窗口下二者的匹配精度。針對粒子均勻移動的理想情況,在相同匹配窗口時,x方向和y方向的匹配精度差別不大。二次多項式曲面擬合法匹配精度隨著匹配窗口的增大而提高,基于梯度亞像素位移算法匹配精度幾乎不隨窗口變化,但兩種亞像素算法的匹配精度均在0.1個像素內(nèi)。

2)匹配窗口為16像素×16像素時,最合適的粒子濃度區(qū)間為0.040～0.360;匹配窗口為32像素×32像素時,最合適的粒子濃度區(qū)間為0.008～0.440;匹配窗口為64像素×64像素時,最合適的粒子濃度區(qū)間為0.002～0.530。匹配窗口越大,最合適的粒子濃度范圍越大。在不同的匹配窗口下,粒子直徑為3個像素時的誤差均較小,為保證粒子跟隨性,建議在PIV實際使用過程中,選擇直徑為3個像素的粒子為宜。

3)針對理想旋渦的粒子圖像,利用基于梯度亞像素位移算法得到的誤差比利用二次多項式曲面擬合法得到的誤差要小,二次多項式曲面擬合法的匹配誤差基本穩(wěn)定在0.2個像素內(nèi),而基于梯度亞像素位移算法得到的誤差基本在0.05個像素以內(nèi)。