劉旭春，李德龍，，劉 津，，劉 楊，，滕長勝，靖常峰

(1.北京建筑工程學院，北京100044;2.現(xiàn)代城市測繪國家測繪地理信息局重點實驗室，北京100044;3.中測新圖(北京)遙感技術有限責任公司，北京100039)

一 、引 言

低空航攝能夠獲取大比例尺高精度的影像，且具有響應時間短和靈活多變的優(yōu)勢，因此常用于應急災害測繪。同時，救災的緊迫性要求遙感影像能夠被快速處理。在快速處理領域，一般采用分布式集群計算機和多核CPU來獲得大數(shù)據(jù)量和復雜計算的加速效果，但其局限性比較大，大集群占用太大空間，多核CPU性能難以達到要求，而可編程圖形處理器(graphic process unit，GPU)同時解決了分布式集群和多核CPU的局限性。尤其是近年來，可編程圖形處理器的并行處理在其他領域得到了很好的應用和發(fā)展，因此研究其在遙感影像快速處理領域的應用非常有實際意義。

可編程圖形處理器利用其可并行運算的CUDA(compute unified device architecture)統(tǒng)一計算設備軟硬件架構體系，可在硬件加速和軟件并行設計方面提高圖像處理的計算速度。本文通過對3種常見影像處理算法的CUDA編程設計實現(xiàn)，得到了在精度不變的情況下，GPU并行比CPU串行算法執(zhí)行速度有顯著提高的結論，大幅度提高了現(xiàn)階段低空遙感影像處理的速度。

二、可編程圖形處理器及編程模型

可編程圖形處理器是由NVIDIA公司開發(fā)的具有高性能數(shù)據(jù)處理特性的集成到顯卡上的專用硬件，它的型號有 TESLA、QUADRO、GTX等。其中，TESLA專用于數(shù)值計算;QUADRO專用于圖像處理和顯示。FREMI 2.0及以上的GPU已經(jīng)基本支持C++和數(shù)值double類型，可以提供更加便利的開發(fā)環(huán)境和更高的精度。GPU硬件正逐漸走向成熟，在硬件架構上主要是采用了統(tǒng)一處理架構和引入了片內(nèi)共享存儲器，支持隨機寫入和線程間通信，因此它可以輕松實現(xiàn)Tflops級計算精度。

CUDA是基于CPU和GPU的異構編程模型，它把CPU作為主機端(host)，GPU作為協(xié)處理器或者設備端(device)，因此在一個系統(tǒng)中可以有一個主機端和若干個設備端，它們協(xié)同合作。在CUDA異構模型中，CPU主要負責邏輯性強的事務處理及與硬盤的數(shù)據(jù)交換，GPU則專注于執(zhí)行高度密集線程化的數(shù)據(jù)并行處理，它們之間的數(shù)據(jù)交換通過內(nèi)存和顯存的數(shù)據(jù)拷貝來完成。PCI-E3.0接口的上下行速率已經(jīng)達到了16 GB/s，GPU的顯存帶寬已經(jīng)達到了140 GB/s。

CUDA編程模型的計算步驟基本分為3步：①待處理數(shù)據(jù)從內(nèi)存拷貝到顯存;②GPU并行數(shù)據(jù)處理;③處理后的數(shù)據(jù)從顯存拷貝回內(nèi)存。GPU多線程處理是通過__global__函數(shù)來驅(qū)動的，稱為內(nèi)核函數(shù)。內(nèi)核函數(shù)的參數(shù)規(guī)定了 Grid、Block與Thread間的層次并行模型。

CUDA編程模型并不一定適用于所有的編程算法，但適合的算法必須滿足它的統(tǒng)一編程模型架構。CUDA編程模型對算法的要求主要有以下幾點：

1)數(shù)據(jù)量很大，計算復雜度高。因為內(nèi)存和顯存的拷貝同樣占用時間，如果數(shù)據(jù)量不大，計算復雜度小，有可能出現(xiàn)CPU比GPU運算速度快的情況。

2)數(shù)據(jù)可分塊。數(shù)據(jù)塊之間具有相互獨立性且計算步驟十分相近。

3)算法本身應盡可能地少用邏輯判斷。因為GPU有很多計算單元，但邏輯單元非常少。

4)GPU編程時一定要注意其數(shù)據(jù)訪存模式。必須實現(xiàn)對全局存儲器的合并訪存，以及對常數(shù)存儲器、紋理存儲器和共享存儲器的合理應用。

三、影像處理算法分析CUDA設計

1.正射糾正算法分析設計

影像正射糾正算法用于改正航攝過程中由航攝軸傾斜、中心投影、大氣折光、地球曲率等因素的影響而產(chǎn)生的影像幾何變形。它的主要流程依次包括建立糾正后影像框架、共線方程計算像點坐標、灰度內(nèi)插賦值。正射糾正算法在CPU端主要的費時操作是共線方程計算像點坐標和灰度插值，這兩個過程有大量的循環(huán)和計算。

1)共線方程可以解算出正射后影像每個點的地面坐標所對應源影像上的每個點像元坐標的映射關系。其公式為

式中，Zuv為對應點的高程，由DEM內(nèi)插求得。

2)灰度內(nèi)插就是把對應的源影像上相關位置的灰度值賦值給正射后的影像，擬采用雙線性內(nèi)插方法。

根據(jù)以上算法分析，每一個像素點都要進行共線方程解算新位置并插值計算像素值。共線方程映射關系具有相對的數(shù)據(jù)獨立性，利用GPU并行處理技術更適合CUDA編程模型。

根據(jù)正射糾正算法分析和CUDA編程模型，整個算法的實現(xiàn)流程如下：①CPU讀取影像和DEM等數(shù)據(jù)，以及糾正后建立的影像，將這些數(shù)據(jù)均存放在內(nèi)存中;②CPU把算法需要的數(shù)據(jù)由內(nèi)存拷貝到顯存，同時為糾正后影像開辟顯存空間，原始影像和DEM數(shù)據(jù)訪問具有不對稱性，可放入紋理存儲器來縮短訪問時間，糾正所需的常用參數(shù)則放入共享存儲器來加速訪問;③在GPU端，利用CUDA編程模型的多線程并行處理技術來實現(xiàn)共線方程，灰度插值在各線程內(nèi)部進行，全局存儲器均采用合并訪存的方式訪問，以避免各項訪問沖突;④在CPU端，把糾正后影像拷貝到內(nèi)存，釋放各種GPU存儲器資源，并將糾正后數(shù)據(jù)寫入硬盤。

2.快速傅里葉變換算法分析

快速傅里葉變換(fast Fourier transform，F(xiàn)FT)的基本思想是利用DFT旋轉(zhuǎn)因子的周期性和對稱性，使長序列的DFT分解為更小長度的DFT，然后利用這些小的DFT計算來代替大的DFT計算，從而達到提高效率的目的。FFT算法在CPU端主要的費時操作是因為像幅大、循環(huán)次數(shù)多、計算較為復雜。

設x(n)為長度為N的有限長序列，則其離散傅里葉正變換(DFT)為

式中，WN=，稱為旋轉(zhuǎn)因子或蝶形因子。

通過兩次調(diào)用一維快速傅里葉變換即可實現(xiàn)二維快速傅里葉變換。二維快速傅里葉變換的輸入相當于一個M行N列的復數(shù)矩陣，對其每行進行一維快速傅里葉變換，再對其每列進行一維快速傅里葉變換，這樣就達到了變換的目的。

根據(jù)以上算法分析，每一個像素都要依次進行行列各一次DFT變換。此算法具有相對的數(shù)據(jù)獨立性，可以利用GPU并行處理技術，適合CUDA編程模型。

本試驗編程調(diào)用了CUFFT函數(shù)庫。CUDA提供了封裝好的CUFFT庫，它提供了與CPU上的FFTW庫相似的接口，能夠讓使用者輕易地挖掘GPU的強大浮點處理能力，不用自己去實現(xiàn)專門的FFT內(nèi)核函數(shù)。使用者通過調(diào)用CUFFT庫的API函數(shù)，即可完成FFT變換。

3.高斯差分算法分析設計

高斯差分是SIFT匹配算法的關鍵步驟，它包括影像空間域卷積和生成DOG影像。影像空間域卷積算法簡單，即用一個窗口模板移動來平滑整幅影像;再把不同卷積模板卷積的影像對應位置相減即得到高斯差分圖像。高斯差分算法在CPU端所需時間主要還是消耗在大量的循環(huán)和復雜的計算中。

本次試驗算法決定使用5幅卷積圖像生成1組4層高斯差分影像。為了簡便算法，設計把二維卷積變換為先行卷積再列卷積，再利用卷積圖像相減得到差分圖像。

根據(jù)以上算法分析，每一個像素點都要先行卷積再列卷積，再對應位置相減。此過程具有天然的數(shù)據(jù)獨立性，非常適合GPU并行CUDA編程模型。

根據(jù)高斯差分算法分析和CUDA編程模型，整個算法的實現(xiàn)流程如下：①CPU讀入影像，生成卷積模板;②CPU把影像和模板數(shù)據(jù)開辟，并拷貝顯存空間，影像放入全局存儲器，采用合并訪存的方式訪問，因為模板具有重復使用性和不改變性，所以放入共享存儲器來提高訪問速度，另外還需在全局存儲器中開辟5幅中間影像空間和5幅高斯差分結果空間;③在GPU端，為了確保列變換之前行變換已完成，采用CUDA中的流控制來依次啟動行列兩個卷積Kernel函數(shù)，利用CUDA編程模型的多線程并行處理技術來拆分多次循環(huán)，卷積在各線程內(nèi)部進行，每個線程負責一個像素點的相減;④在CPU端把4幅高斯差分影像拷貝到內(nèi)存，釋放各種GPU存儲器資源，并將差分影像寫入硬盤。

四、CUDA試驗及精度時間分析

1.試驗環(huán)境

為了更有效地說明圖形處理器在遙感影像快速處理領域的作用，本文所有編程試驗均在統(tǒng)一的電腦配置環(huán)境下進行，除了添加一個GPU Quadro 4000顯卡外，其余配置都和經(jīng)典編程環(huán)境一致。

算法實現(xiàn)的主要硬件平臺：CPU為Intel(R)Xeon(R)CPU W3550@3.07 GHz4核8線程;GPU為NVIDIA Quadro 4000。

算法依賴的主要軟件平臺：Windows XP X86、Visual Studio 2008、NVIDIA GPU Computing Toolkit(CUDA 4.1)、FFTW和 CUFFT函數(shù)庫。

為了更好地體現(xiàn)圖形處理器的性能，不僅要對多幅不同像幅的原始遙感影像進行處理，還要分析算法在CPU和GPU上運算數(shù)據(jù)精度和時間的差異。

2.CPU與GPU數(shù)據(jù)精度分析

由于CPU和GPU由不同的處理單元組成，對單浮點和雙浮點類型的支持可能不同，因此有可能會造成除法保留的小數(shù)位不同，最終造成影像灰度級別的差異。同時，由于算法必須為適合CUDA統(tǒng)一編程模型而作適當?shù)男薷?，同樣可能存在差異，因此必須對比它們的?shù)據(jù)精度。在此，僅以單波段正射糾正為例，將基于CPU和GPU獲得的影像進行數(shù)據(jù)精度分析，統(tǒng)計兩者的差異。

圖1與圖3分別表示CPU正射糾正后部分影像圖和部分像素值圖(保留小數(shù)點后3位)，圖2與圖4則表示GPU的相同內(nèi)容。圖1和圖2對比，影像圖在視覺上無差別;圖3和圖4對比，相對應部分像素值內(nèi)容完全相同。為了確保結論的正確性，令CPU與GPU正射糾正后的兩幅影像相對應像素值相減，再對結果進行統(tǒng)計，發(fā)現(xiàn)所有對應像素位置的差值均為零。綜上所述，說明GPU在數(shù)據(jù)精度上與CPU處理精度相同，達到了遙感影像快速處理對數(shù)據(jù)精度的要求。

圖1 CPU正射糾正后部分影像圖

圖2 GPU正射糾正后部分影像圖

圖3 CPU影像部分像素值

圖4 GPU影像部分像素值

3.CPU與GPU時間分析

為了更加充分地說明GPU加速遙感影像處理的能力，采用像幅大小分別為1024像素×1024像素、2048像素×2048像素、4096像素×4096像素的遙感影像進行時間統(tǒng)計對比分析。測試結果見表1。

由表1計算可見，正射糾正的最高加速比為1947/377=5.16;FFT的最高加速比為2453/233=10.53;高斯差分的最高加速比為3057/330=9.26。由圖5分析可得，GPU處理遙感影像的加速比隨著像幅的增大而增加，這說明在不超過GPU存儲單元容量的條件下，遙感影像數(shù)據(jù)量越大越能體現(xiàn)GPU的加速效果，也能得到更大的加速比。僅從表1處理時間上分析，利用GPU并行處理技術可以在4096像素×4096像素像幅的情況下獲得比CPU串行高達一個量級的加速比，向影像快速處理邁進了重要一步。

表1 測試結果ms

圖5 GPU加速比

五、結束語

基于以上的試驗結果分析，GPU能夠在影像正射糾正、快速傅里葉變換和高斯差分方面得到明顯的加速，并能夠在數(shù)據(jù)精度上與CPU保持一致的情況下，大幅提高遙感影像的處理速度。遙感影像處理是針對像素的操作，數(shù)據(jù)處理具有以下特點：數(shù)據(jù)量大、計算復雜、處理方法固定、處理方法具有內(nèi)在的并行特性等。遙感影像處理的這些特點決定了它適合于GPU的并行處理架構和CUDA編程模型。筆者認為利用GPU技術能夠?qū)崿F(xiàn)遙感圖像快速處理領域的巨大飛躍。

［1］張舒，楮艷利.GPU高性能運算之CUDA［M］.北京：中國水利水電出版社，2009：10.

［2］張劍清，潘勵，王樹根.攝影測量學［M］.武漢：武漢大學出版社，2003：158-162.

［3］周海芳，遙感圖像并行處理算法的研究與應用［D］.長沙：國防科技大學，2003.

［4］劉勇，使用GPU加速通用科學計算－CUDA技術解析［J］.科技信息，2008(24)：394-395.

［5］楊靖宇，張永生，張宏蘭，等.基于可編程圖形硬件的遙感影像并行處理研究［J］.測繪工程，2008，17(3)：21-22.

［6］楊志義，朱婭婷，薄勇，基于統(tǒng)一計算設備架構技術的并行圖像處理研究［J］.計算機測量與控制，2009，17(4)：734-737.

［7］SINHA SN，F(xiàn)RAHM JM，POLLEFEYSM，et al.Feature Tracking and Matching in Video Using Programmable Graphics Handware［J］.Machine Vision and Application，2007，22(1)：207-217.

［8］NVIDIA.NVIDIA CUDA Computer Unified Device Architecture Programming Guide［M］.version 2.0.NVIDIA：［s.n.］，2008.

［9］GARCIA V，DEBREUVE E，BARLAUD M.Fast KNearest Neighbor Search using GPU［C］∥IEEE Confidence on Computer Version and Patter Recognition Workshops.Alaska：［s.n.］，2008：1-6.

［10］VASILIADIS G，ANTONATOS S，POLYCHRONAKIS M，et al.Gnort：High Performance Network Intrusion Detection Using Graphics Processors［C］∥Recent Advance in Intrusion Detection(RAID).Heidelberg：Springer-Verlay，2008：116-134.