王移風 曹衍龍 徐 朋 金 鷺 楊將新

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

0 引言

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

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

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

1 系統(tǒng)結構

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

圖1 系統(tǒng)結構

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

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

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

可得點（x，y）的絕對光程差h（x，y）為

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

本文采用快速傅里葉變換法［11-12］來對光強分布函數進行分析，此時式（1）可改寫為

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

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

式（7）中的虛數部分即為每個波數所對應的相位值。此時的相位值是以2π為周期分布的，并不能直接應用，需要先進行相位的展開運算，通過展開的相位值可以求解該點的絕對光程差。

2 基于GPU的數據處理算法

鑒于上述算法在處理干涉圖像時，是對干涉圖上的各像素點光強進行分析，而每個像素點之間的處理是相互獨立的，故可以采用并行處理的方式對這一過程進行加速。

2.1 統(tǒng)一計算設備架構CUDA［12］

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

2.2 基于CUDA的GPU并行算法的設計

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

圖2 傅里葉變換分析法的GPU并行實現流程圖

（1）對光強分布函數進行快速傅里葉變換。首先將原始干涉信號進行快速傅里葉變換。

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

（3）傅里葉反變換。濾去a（σ）之后，利用CUDA做傅里葉反變換，將信號從頻域轉換到時域上，此時存在歸一化問題，即傅里葉反變換后的值應除以上步FFT的點數才是其真實值。

（4）求取相位信息。將過濾a（σ）后的信號做對數運算，得到結果的虛數部分即為相位值。這一過程對于每一個波數σ都是獨立的，故可以選擇大尺度的線程塊維度。

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

圖3 過濾前后的幅值分布

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

式中，k為相位周期數。

圖4 相位展開原理

（7）求取絕對光程差。展開相位后，得到每個波數所對應的相位值，由最小二乘法擬合出相位分布，根據式（4）求取該像素點對應工件表面的絕對光程差。

3 實驗結果分析

3.1 結果比較

基于GPU的相位計算并行處理算法與CPU運行下的CPU串行處理算法進行實驗比較。

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

圖5c表明，CPU和GPU運算結果相對誤差最大僅為1.5×10－6，可以忽略不計。

3.2 加速對比

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

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

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

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

4 結束語

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

［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，Fexon L，Feruglio 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.