劉旭春，李德龍，，劉 津，，劉 楊，，滕長(zhǎng)勝，靖常峰

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

一 、引 言

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

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

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

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

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

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

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

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

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

3)算法本身應(yīng)盡可能地少用邏輯判斷。因?yàn)镚PU有很多計(jì)算單元，但邏輯單元非常少。

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

三、影像處理算法分析CUDA設(shè)計(jì)

1.正射糾正算法分析設(shè)計(jì)

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

1)共線方程可以解算出正射后影像每個(gè)點(diǎn)的地面坐標(biāo)所對(duì)應(yīng)源影像上的每個(gè)點(diǎn)像元坐標(biāo)的映射關(guān)系。其公式為

式中，Zuv為對(duì)應(yīng)點(diǎn)的高程，由DEM內(nèi)插求得。

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

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

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

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

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

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

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

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

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

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

3.高斯差分算法分析設(shè)計(jì)

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

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

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

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

四、CUDA試驗(yàn)及精度時(shí)間分析

1.試驗(yàn)環(huán)境

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

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

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

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

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

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

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

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

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

圖3 CPU影像部分像素值

圖4 GPU影像部分像素值

3.CPU與GPU時(shí)間分析

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

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

表1 測(cè)試結(jié)果ms

圖5 GPU加速比

五、結(jié)束語(yǔ)

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

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

［2］張劍清，潘勵(lì)，王樹(shù)根.攝影測(cè)量學(xué)［M］.武漢：武漢大學(xué)出版社，2003：158-162.

［3］周海芳，遙感圖像并行處理算法的研究與應(yīng)用［D］.長(zhǎng)沙：國(guó)防科技大學(xué)，2003.

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

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

［6］楊志義，朱婭婷，薄勇，基于統(tǒng)一計(jì)算設(shè)備架構(gòu)技術(shù)的并行圖像處理研究［J］.計(jì)算機(jī)測(cè)量與控制，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.