王移風(fēng) 曹衍龍 徐 朋 金 鷺 楊將新

1.浙江科技學(xué)院，杭州，310023 2.浙江理工大學(xué)，杭州，310018 3.浙江大學(xué)，杭州，310027

0 引言

隨著機(jī)械加工精度的不斷提高，與表面形貌有關(guān)的應(yīng)用也越來越廣泛，表面形貌測量技術(shù)在精密加工行業(yè)和各類研究領(lǐng)域中的地位顯得愈發(fā)重要。白光譜線掃描干涉法是表面形貌檢測的有效方法，適合于大范圍、高精度的測量，但由于圖像數(shù)據(jù)量大，數(shù)據(jù)處理時(shí)間過長，限制了在實(shí)時(shí)測量中的應(yīng)用。如何加速數(shù)據(jù)處理的過程，縮短數(shù)據(jù)處理時(shí)間，提升系統(tǒng)整體的實(shí)時(shí)性已成為關(guān)鍵問題。

圖形處理器GPU可以在較低的成本下，以較快的速度完成科學(xué)計(jì)算中的并行處理問題，并具有十分明顯的加速效果，近年來圖形處理器已在諸如光學(xué)測量［1］、物理模擬［2－4］、科學(xué)計(jì)算［5－6］、圖像處理［7－10］等領(lǐng)域得到成功應(yīng)用。

本文介紹一種利用GPU進(jìn)行數(shù)據(jù)處理的實(shí)時(shí)表面形貌測量方法，利用NVIDIA Geforce GTS250圖形卡（16個(gè)流多處理器，1100MHz）進(jìn)行數(shù)據(jù)處理，對比分析了CPU運(yùn)行下的串行處理算法的運(yùn)行時(shí)間和精度，實(shí)驗(yàn)表明基于GPU的數(shù)據(jù)處理可以得到很好的效果。

1 系統(tǒng)結(jié)構(gòu)

所搭建的表面形貌測量系統(tǒng)主要由白光譜線掃描干涉模塊、聲光濾光器控制模塊、圖像采集模塊和數(shù)據(jù)處理模塊組成。其系統(tǒng)結(jié)構(gòu)如圖1所示。白光源經(jīng)過透鏡發(fā)出平行光，經(jīng)過濾光裝置后成為波長呈線性變化的單色光。這束單色光進(jìn)入類似邁克爾遜干涉儀的系統(tǒng)之中，經(jīng)過分光鏡BS后分為兩束，一束打在參考反射鏡上，另一束打在被測工件的表面。兩束光反射后，重新在分光鏡BS處匯聚并產(chǎn)生干涉條紋。干涉圖像由CCD采集并存放到計(jì)算機(jī)中進(jìn)行處理。

圖1 系統(tǒng)結(jié)構(gòu)

干涉圖像上單個(gè)像素點(diǎn)的灰度值表征了某一波長下該點(diǎn)發(fā)生干涉時(shí)其干涉條紋的光強(qiáng)信息。由CCD上某像素點(diǎn) （x，y）探測到的光強(qiáng)可表示為

式中，I1（σ）、I2（σ）分別為由被測工件與參考反射鏡反射回來的光線光強(qiáng)；σ為波數(shù)；φ（x，y，σ）為干涉信號(hào)的相位；h（x，y）為該點(diǎn)處兩束反射光的絕對光程差。

由波數(shù)改變而引起的相位改變?yōu)?/p>

可得點(diǎn)（x，y）的絕對光程差h（x，y）為

濾光裝置可以使得波數(shù)呈線性變化，同時(shí)CCD的采集過程也是線性的，因此通過求取波數(shù)變化Δσ時(shí)的相位變化Δφ（x，y，k），就能求得該點(diǎn)的絕對光程差。通過分析整張干涉圖像上所有像素點(diǎn)的相位變化得到絕對光程差，以此能描繪出該被測區(qū)域表面的形貌特征，以進(jìn)行下一步的分析和評定，由此可見精確的相位是獲得絕對光程差的關(guān)鍵。

本文采用快速傅里葉變換法［11-12］來對光強(qiáng)分布函數(shù)進(jìn)行分析，此時(shí)式（1）可改寫為

a（σ）和b（σ）相對于光程差h（x，y）來說是緩變元，其頻率遠(yuǎn)低于后者，故可以在頻域上對低頻部分進(jìn)行濾波，將a（σ）濾去。對式（5）做傅里葉變換：

將a（σ）過濾之后做反傅里葉變換，并通過對數(shù)運(yùn)算分離相位信息：

式（7）中的虛數(shù)部分即為每個(gè)波數(shù)所對應(yīng)的相位值。此時(shí)的相位值是以2π為周期分布的，并不能直接應(yīng)用，需要先進(jìn)行相位的展開運(yùn)算，通過展開的相位值可以求解該點(diǎn)的絕對光程差。

2 基于GPU的數(shù)據(jù)處理算法

鑒于上述算法在處理干涉圖像時(shí)，是對干涉圖上的各像素點(diǎn)光強(qiáng)進(jìn)行分析，而每個(gè)像素點(diǎn)之間的處理是相互獨(dú)立的，故可以采用并行處理的方式對這一過程進(jìn)行加速。

2.1 統(tǒng)一計(jì)算設(shè)備架構(gòu)CUDA［12］

在CUDA架構(gòu)中，CPU作為主機(jī)執(zhí)行一些復(fù)雜的邏輯處理和事務(wù)處理，同時(shí)作為協(xié)處理器的設(shè)備端GPU負(fù)責(zé)執(zhí)行高度線程化的并行處理任務(wù)。CUDA程序就是由主機(jī)端的串行執(zhí)行部分加上設(shè)備端并行執(zhí)行的一系列內(nèi)核函數(shù)組成的。內(nèi)核函數(shù)通過一個(gè)線程網(wǎng)格（Grid）來執(zhí)行。執(zhí)行內(nèi)核函數(shù)的線程網(wǎng)格包含著兩層并行的結(jié)構(gòu)，即線程塊（Block）和線程（Thread），線程塊之間以及每個(gè)線程塊中的所有線程之間都是相互獨(dú)立并可以并行執(zhí)行程序的。本文的并行處理算法就是通過構(gòu)建一系列內(nèi)核函數(shù)實(shí)現(xiàn)的。

2.2 基于CUDA的GPU并行算法的設(shè)計(jì)

傅里葉變換分析法的GPU并行實(shí)現(xiàn)流程如圖2所示。首先，在內(nèi)存中完成初始化；其次，拷貝到顯存運(yùn)用FFT－IFFT方法濾除多余的光強(qiáng)信號(hào)并提取光強(qiáng)分布函數(shù)的相位信息；最后，通過求解相位跳變并展開干涉信號(hào)相位求取絕對光程差。

圖2 傅里葉變換分析法的GPU并行實(shí)現(xiàn)流程圖

（1）對光強(qiáng)分布函數(shù)進(jìn)行快速傅里葉變換。首先將原始干涉信號(hào)進(jìn)行快速傅里葉變換。

（2）帶通濾波。設(shè)計(jì)內(nèi)核函數(shù)進(jìn)行帶通濾波將式（5）中的變元a（σ）過濾掉。過濾前后的頻率分布如圖3所示。

（3）傅里葉反變換。濾去a（σ）之后，利用CUDA做傅里葉反變換，將信號(hào)從頻域轉(zhuǎn)換到時(shí)域上，此時(shí)存在歸一化問題，即傅里葉反變換后的值應(yīng)除以上步FFT的點(diǎn)數(shù)才是其真實(shí)值。

（4）求取相位信息。將過濾a（σ）后的信號(hào)做對數(shù)運(yùn)算，得到結(jié)果的虛數(shù)部分即為相位值。這一過程對于每一個(gè)波數(shù)σ都是獨(dú)立的，故可以選擇大尺度的線程塊維度。

（5）求取相位跳變。由FFT－IFFT法求得的相位呈周期性分布，需對其進(jìn)行相位展開。在展開之前先求取其相位跳變，分析其相位分布情況。

圖3 過濾前后的幅值分布

（6）相位展開。相位展開原理如圖4所示。其表達(dá)形式為

式中，k為相位周期數(shù)。

圖4 相位展開原理

（7）求取絕對光程差。展開相位后，得到每個(gè)波數(shù)所對應(yīng)的相位值，由最小二乘法擬合出相位分布，根據(jù)式（4）求取該像素點(diǎn)對應(yīng)工件表面的絕對光程差。

3 實(shí)驗(yàn)結(jié)果分析

3.1 結(jié)果比較

基于GPU的相位計(jì)算并行處理算法與CPU運(yùn)行下的CPU串行處理算法進(jìn)行實(shí)驗(yàn)比較。

圖5a、圖5b所示為150像素×150像素區(qū)域CPU和GPU的運(yùn)算結(jié)果；圖5c為該塊絕對光程差數(shù)據(jù)的相對誤差，其定義為：相對誤差＝（GPU結(jié)果－CPU結(jié)果）／CPU結(jié)果。

圖5c表明，CPU和GPU運(yùn)算結(jié)果相對誤差最大僅為1.5×10－6，可以忽略不計(jì)。

3.2 加速對比

設(shè)計(jì)開發(fā)了基于GPU的相位計(jì)算并行處理算法，并與CPU運(yùn)行下的CPU串行處理算法的運(yùn)行時(shí)間進(jìn)行實(shí)驗(yàn)比較，見表1。為對比不同數(shù)據(jù)量下的兩種算法的運(yùn)算效率，分別對多個(gè)不同大小范圍的圖像區(qū)域進(jìn)行處理。處理分為相位求取階段和相位展開求光程差階段。

表1 傅里葉分析法GPU及CPU處理速度及加速比

由表1可以得出：GPU加速比隨著數(shù)據(jù)處理數(shù)量的增多而增大。在數(shù)據(jù)量小的情況下內(nèi)存和顯存之間數(shù)據(jù)搬運(yùn)的時(shí)間沒有被GPU運(yùn)算時(shí)間掩蓋，使用CUDA進(jìn)行小計(jì)算量的運(yùn)算是不劃算的，只有當(dāng)并行計(jì)算在整個(gè)應(yīng)用中所占的比例較大時(shí)，其加速性能才能很好地得到體現(xiàn)。

相位展開階段的加速比大于相位求取階段的加速比，這是由于相位展開階段全部采用速度極快的共享存儲(chǔ)器進(jìn)行并行運(yùn)算，并且在這個(gè)過程中并行算法避免了程序出現(xiàn)共享存儲(chǔ)器區(qū)沖突問題。

4 結(jié)束語

本文基于CUDA提出了表面形貌測量的數(shù)據(jù)處理算法，經(jīng)實(shí)驗(yàn)表明運(yùn)算速度比傳統(tǒng)的算法大幅提高，其運(yùn)算精度也可以得到保證?；贕PU的數(shù)據(jù)處理算法對表面形貌實(shí)時(shí)測量有重要意義。

［1］Zhang S.GPU－assisted High－resolution，Realtime 3－D Shape Measurement［J］.Optics Express，2006，14（20）：9120－9130.

［2］Hu L，Xuan L，Li D，et al.Real－time Liquid－crystal Atmosphere Turbulence Simulator with Graphic Processing Unit［J］.Optics Express，2009，17（9）：7259－7269.

［3］Ren N，Liang J，Qu X，et al.GPU－based Monte Carlo Simulation for Light Propagation in Complex Heterogeneous Tissues［J］.Optics Express，2010，18（7）：6811－6824.

［4］Fang Q，Boas D A.Monte Carlo Simulation of Photon Migration in 3DTurbid Media Accelerated by Graphics Processing Units［J］.Optics Express，2009，17（22）：20178－20192.

［5］Kang H，YamaguchiT，Yoshikawa H，et al.Acceleration Method of Computing A Compensated Phase－added Stereogram on A Graphic Processing Unit［J］.Applied Optics，2008，47（31）：5784－5790.

［6］Bin H，Tarek T.Acceleration of Spiking Neural Network Based Pattern Recognition on NVIDIA Graphics Processors［J］.Applied Optics，2010，49（10）：B83－B91.

［7］Pan Y，Xu X，Solanki S，et al.Fast CGH Computation Using S－LUT on GPU［J］.Optics Express，2009，17（21）：18543－18556.

［8］Kang H，Yaras F，Onural L.Graphics Processing U－nit Accelerated Computation of Digital Holograms［J］.Appled Optics，2009，48（34）：H137－H144.

［9］Vinegoni C，F(xiàn)exon L，F(xiàn)eruglio P F，et al.High Throughput Transmission Optical Projection Tomography Using Low Cost Graphics Processing Unit［J］.Optics Express，2009，17（25）：22320－22333.

［10］Shimobaba T，Sato Y，Miura J，et al.Real－time Digital Holographic Microscopy Using the Graphic Processing Unit［J］.Optics Express，2008，16（16）：11776－11782.

［11］Sandoz P，Tribillon G，Perrin H.High－resolution Profilometry by Using Phase Calculation Algorithms for Spectroscopic Analysis of White－light Interferograms［J］.Journal of Modern Optics，1996，43（4）：701－708.

［12］Kuwamura S，Yamaguchi I.Wavelength Scanning Profilometry for Real－time Surface Shape Measurement［J］.Applied Optics，1997，36：4473－4482.