一種改進的基于FFT的PIV互相關算法

鮑曉利， 李木國

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

0 引 言

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

1 粒子圖像測速原理

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

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

2 基于FFT的互相關算法

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

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

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

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

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

3 基于FFT的互相關算法的改進

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

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

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

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

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

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

4 實驗結果與討論

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

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

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

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

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

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

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

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

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

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

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

5 結 論

本文利用頻域抽取原理，在基于FFT的互相關算法的基礎上，結合重疊窗口互相關運算的特點，對基于FFT的互相關算法進行了改進，有效地減少了互相關運算中重復FFT運算量.實驗結果表明，本文所提出的改進互相關算法與基于FFT的互相關算法相比，在不降低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］.水動力學研究與進展，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技術及其在流場測量中的應用［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］胡廣書.數字信號處理（理論、算法與實現）［M］.2版.北京：清華大學出版社，2003：176－177