一種改進的基于FFT的PIV互相關(guān)算法

鮑曉利， 李木國

（大連理工大學(xué)海岸和近海工程國家重點實驗室，遼寧大連 116024）

0 引 言

近年來，粒子圖像測速（particle image velocimetry，PIV）技術(shù)得到了極大的發(fā)展，已被廣泛應(yīng)用于多個科研和工業(yè)領(lǐng)域.在2D－PIV技術(shù)中，基于相關(guān)窗口的互相關(guān)算法已被普遍采用［1］.隨著對3D－PIV技術(shù)研究的不斷深入，基于相關(guān)窗口的互相關(guān)算法也被應(yīng)用其中［2、3］.在現(xiàn)有的實時PIV系統(tǒng)設(shè)計中，也多采用基于相關(guān)窗口的互相關(guān)算法［4］.同時，基本相關(guān)窗口的互相關(guān)算法也被用于相關(guān)領(lǐng)域，如三維運動檢測［5］.當(dāng)前對互相關(guān)算法的研究主要集中在提高算法精度方面，在提高算法運算效率方面的研究沒有得到足夠的重視，所以對互相關(guān)算法的運算效率進行研究，具有重要意義.到目前為止，在所提出的互相關(guān)算法中，基于快速傅里葉變換（FFT）的互相關(guān)算法由于其不受相關(guān)窗口大小限制，簡化了互相關(guān)計算，并極大地提高了互相關(guān)運算效率，已成為互相關(guān)算法中最為典型的算法［6、7］.在粒子圖像測速中，為提高粒子速度矢量分辨率，相鄰相關(guān)窗口間須存在重疊，致使互相關(guān)運算中存在重復(fù)運算，降低了算法的運算效率［8］.本文在基于FFT的互相關(guān)算法的基礎(chǔ)上，根據(jù)頻域抽取原理，對互相關(guān)算法中的二維FFT進行改進，以減少互相關(guān)算法中的重復(fù)FFT運算.

1 粒子圖像測速原理

PIV技術(shù)是在流動顯示的基礎(chǔ)上結(jié)合光學(xué)圖像處理技術(shù)對待測量的流場進行非接觸的高精度測量［9］.測速系統(tǒng)對加入示蹤粒子的流場進行快速拍照，通過測量流場中示蹤粒子在已知時間間隔內(nèi)的位移實現(xiàn)對流體運動速度的測量.定義x（t）、y（t）為某個粒子在t時刻的位置，經(jīng)過時間Δt粒子的位置為x（t＋Δt）、y（t＋Δt），那么該粒子在x軸、y軸方向的速度分量v x、v y如下式：

實際應(yīng)用中，粒子在t＋Δt時刻的位置是未知的.對于粒子在t＋Δt時刻的位置，通常采用基于相關(guān)窗口的互相關(guān)方法來確定.

2 基于FFT的互相關(guān)算法

圖1 基于FFT的互相關(guān)算法流程圖Fig.1 FFT－based cross－correlation algorithm diagram

針對互相關(guān)算法運算量巨大的問題提出了較多算法，基于FFT的互相關(guān)算法由于具有極高的運算效率，至今仍為PIV互相關(guān)算法中的基本算法［10］.圖1為基于FFT的互相關(guān)算法流程圖.其中f（m，n）和g（m，n）分別表示當(dāng)前窗口和詢問窗口.t時刻粒子存在于（x，y）位置處窗口（即當(dāng)前窗口），t＋Δt時刻粒子所在窗口未知.t＋Δt時刻粒子所在窗口位置（x＋ξ，y＋η）可利用當(dāng)前窗口與t＋Δt時刻的粒子存在目標區(qū)域進行互相關(guān)運算獲得.最大互相關(guān)峰值位置處窗口，即t＋Δt時刻粒子所在窗口.當(dāng)前窗口與詢問窗口的互相關(guān)函數(shù)如式（2）［11］.其中，當(dāng)前窗口大小為M0×N0，

I1和I2分別表示在t時刻和t＋Δt時刻粒子圖像的像素灰度值.

在基于FFT的互相關(guān)算法中，乘法運算量在互相關(guān)運算量中占絕對比重，所以對互相關(guān)運算量分析可簡化為對其乘法運算量分析，具有等效的分析結(jié)果［11］.基于FFT的互相關(guān)算法的運算量計算表達式為

其中N×N為詢問窗口大小，N滿足2k（k為正整數(shù)），否則詢問窗口須補零.

3 基于FFT的互相關(guān)算法的改進

一般情況下，相關(guān)窗口的重疊率為50%時，粒子速度矢量分辨率與算法運算效率之間達到最佳均衡［8］.頻域抽?。―IF）基二維FFT算法中，信號的FFT由被分成上、下兩部分的信號的FFT構(gòu)成.據(jù)此原理，在窗口重疊率為50%時，對當(dāng)前窗口和詢問窗口進行對半分割，重疊相鄰3個窗口的二維FFT中重疊部分窗口的一個維度的FFT運算值，可由其相鄰窗口的重疊部分的FFT運算值經(jīng)頻移疊加后獲得.據(jù)此，對基于FFT的互相關(guān)算法進行改進，以提高互相關(guān)算法運算效率.在詢問窗口大小為M×N（M、N為2的整數(shù)次冪，一般情況下M等于N）時，重疊窗口的FFT表達式F2（u，v）推導(dǎo)過程如式（4）～（6）［12］所示：

圖2 改進的二維FFT計算流程Fig.2 Computational diagram of improved 2D－FFT

改進后互相關(guān)算法運算量仍采用乘法運算量來度量，其計算表達式如下式：

式（7）與式（3）相比較，改進后算法運算效率提高表達式如式（8）所示.改進后算法較未改進算法運算效率提高如圖3（圖中橫坐標L表示詢問窗口寬度，縱坐標σ表示算法運算效率提高）.

圖3 改進后算法運算效率提高Fig.3 Computation efficiency increase of improved algorithm

從圖3中可以看出，當(dāng)詢問窗口大小為16× 16、32×32、64×64、128×128、256×256、512× 512、1 024×1 024時，改進后算法較未改進算法在運算效率方面平均提高了13.86%，并隨著詢問窗口變大算法運算效率提高會減小.當(dāng)詢問窗口大小為16×16時運算效率提高最大，為15.38%.當(dāng)詢問窗口大小為256×256、512× 512、1 024×1 024時，其運算效率提高變化不明顯.

4 實驗結(jié)果與討論

為對改進前與改進后算法進行對比驗證及分析，利用CCD相機連續(xù)采集了30幀粒子圖像，連續(xù)兩幀粒子圖像時間間隔為0.043 s，圖像大小為1 600×1 200.隨機選取連續(xù)兩幀經(jīng)預(yù)處理后的粒子圖像（為清晰顯示出圖像中粒子，僅給出部分放大圖像），如圖4所示.

當(dāng)詢問窗口大小為16×16或1 024×1 024時，由于詢問窗口所包含粒子信息過少或分析點數(shù)過少，粒子圖像測速處理結(jié)果無任何實際意義.對于這兩種情況本文均不做分析.

設(shè)置窗口重疊率為50%，在詢問窗口大小分別為32×32、64×64、128×128、256×256、512× 512時，采用基于改進前后兩種互相關(guān)算法的PIV算法計算粒子速度矢量.處理程序采用Matlab作為開發(fā)平臺編寫完成.限于篇幅僅給出當(dāng)詢問窗口大小為128×128時隨機選取的粒子速度矢量對比圖，如圖5所示（圖中橫坐標和縱坐標的單位為像素）.

圖4 連續(xù)兩幀粒子圖像Fig.4 A pair of particle images

圖5 粒子速度矢量對比圖Fig.5 Comparison of two particle velocity vector′s images

對比觀察圖5（a）和（b）兩粒子速度矢量圖，發(fā)現(xiàn)其無任何差異.為定量分析改進前后兩算法在所得速度矢量上的差異，對不同詢問窗口大小下得到的粒子速度矢量在x軸和y軸方向上的平均差異進行對比分析，分析結(jié)果如表1所示.

表1 粒子速度矢量對比Tab.1 Comparison results of particle velocity vectors

由表1可以得出，利用改進前后兩互相關(guān)算法計算所得速度矢量在x軸和y軸方向上僅存在微小差異.

改進前后兩種互相關(guān)算法平均消耗時間對比如表2所示.

表2 算法消耗時間對比Tab.2 Comparison results of the processing time

由表2可以得出，改進后算法較未改進算法在運算效率方面平均提高了12.25%，與理論提高效率平均值十分接近.當(dāng)詢問窗口大小為32× 32時改進算法效率提高最高，為14.37%.當(dāng)詢問窗口大小為128×128、256×256、512×512時，隨著詢問窗口變大運算效率提高變化不明顯.

5 結(jié) 論

本文利用頻域抽取原理，在基于FFT的互相關(guān)算法的基礎(chǔ)上，結(jié)合重疊窗口互相關(guān)運算的特點，對基于FFT的互相關(guān)算法進行了改進，有效地減少了互相關(guān)運算中重復(fù)FFT運算量.實驗結(jié)果表明，本文所提出的改進互相關(guān)算法與基于FFT的互相關(guān)算法相比，在不降低PIV處理效果的前提下，運算效率平均提高了約12.25%.

［1］ADIRAN R J.Twenty years of particle image velocimetry［J］.Experiments in Fluids，2005，39（2）：159－169

［2］ZHANG Wei，HAIN R，KHLER C J.Scanning PIV investigation of the laminar separation bubble on a SD7003 airfoil［J］.Experiments in Fluids，2008，45（4）：725－743

［3］WIENEKE B.Volume self－calibration for 3D particle image velocimetry［J］.Experiments in Fluids，2008，45（4）：549－556

［4］MUNOZ J M I，DELLAVALE D，SONNAILLON M O，etal.Real－time particle image velocimetry based on FPGA technology［C］／／5thSouthern Conference on Programmable Logic.Brazil：IEEE，2009

［5］ALVAREZ L，CASTANO C A，CARCIA M，etal.3D motion estimation using a combination of correlation and variational methods for PIV［J］.Computer Aided Systems Theory，2007（4739）：612－620

［6］ADRIC E，VLACHOS P P.Assessment of advanced windowing techniques for digital particle image velocimetry［J］.Measurement Science and Technology，2009，20（7）：1－9

［7］王燦星，林建忠，山本富士夫.二維PIV圖像處理算法［J］.水動力學(xué)研究與進展，2001，16（4）：399－404

［8］ROTH G I，KATZ J.Five techniques for increasing the speed and accuracy of PIV interrogation［J］.Measurement Science and Technology，2001，12（3）：238－245

［9］高殿榮，王益群，申功炘.DPIV技術(shù)及其在流場測量中的應(yīng)用［J］.液壓氣動與密封，2001，21（5）：30－33

［10］STANISLAS M，OKAMOTO K，KHLER C J.Main results of the second international PIV challenge［J］.Experiments in Fluids，2005，39（2）：170－191

［11］TOSHIHITO F，KENJI F，TSUTOMU M.A realtime visualization system for PIV［J］.Field－Programmable Logic and Applications，2003，2778（1）：437－447

［12］胡廣書.數(shù)字信號處理（理論、算法與實現(xiàn)）［M］.2版.北京：清華大學(xué)出版社，2003：176－177