孫業(yè)超 朱曉波 郝雪濤

(中國資源衛(wèi)星應(yīng)用中心，北京 100094)

1 引言

光學(xué)成像系統(tǒng)調(diào)制傳遞函數(shù)(Modulation Transfer Function，MTF)是表征光波經(jīng)成像系統(tǒng)后在頻率域中調(diào)制衰退的函數(shù)，是光學(xué)成像系統(tǒng)性能的一個(gè)重要評(píng)價(jià)指標(biāo)。由于遙感圖像的光學(xué)信息獲取與傳輸過程要經(jīng)過目標(biāo)、大氣、光學(xué)系統(tǒng)、CCD等一系列環(huán)節(jié)，各個(gè)環(huán)節(jié)均可能對(duì)圖像產(chǎn)生退化作用，引起圖像質(zhì)量(Quality)的下降[1]。MTF值反映了成像在不同空間頻率的衰減情況，即頻率越高，信號(hào)衰減的越多，MTF值就越小，所以圖像在采樣數(shù)字化的過程中損失的主要是高頻信息，它在圖像上表現(xiàn)為影像的邊緣、細(xì)節(jié)等的模糊。為了改善成像的清晰度和成像質(zhì)量，需要在地面處理系統(tǒng)中對(duì)數(shù)據(jù)進(jìn)行調(diào)制傳遞函數(shù)補(bǔ)償(Modulation Transfer Function Compensation，MTFC)，恢復(fù)影像的高頻信息，讓影像的邊緣更為清晰,細(xì)節(jié)更為豐富。

隨著遙感對(duì)地觀測技術(shù)的不斷發(fā)展，影像空間分辨率和量化分辨率的提高，使得下傳影像數(shù)據(jù)量不斷增大，目前的MTFC處理研究都是以中央處理器為核心，對(duì)于高空間分辨率的衛(wèi)星影像處理，由于計(jì)算量大造成運(yùn)行速度慢、耗時(shí)長，成為地面系統(tǒng)實(shí)時(shí)處理的瓶頸。

GPU是新一代高性價(jià)比的計(jì)算技術(shù)，近年來發(fā)展迅速。圖形處理器在并行數(shù)據(jù)運(yùn)算上具有強(qiáng)大的運(yùn)算功能以及相對(duì)較高的并行運(yùn)算速度，具有單指令流多線程(Single Instruction Multiple Thread, SIMT)的并行處理特性[2]，在解決計(jì)算密集型問題時(shí)具有較高的性價(jià)比[3]。2006年11月，英偉達(dá)公司推出了計(jì)算統(tǒng)一設(shè)備架構(gòu)(Compute Unified Device Architecture, CUDA)，這是一種新的并行編程模型和指令集架構(gòu)的通用計(jì)算架構(gòu)，能夠利用英偉達(dá) GPU的并行計(jì)算引擎解決復(fù)雜計(jì)算任務(wù)，目前已經(jīng)在醫(yī)學(xué)成像與分割、大氣輻射傳輸計(jì)算、圖像編碼以及光纖通信等領(lǐng)域中得到應(yīng)用，顯著地提高了傳統(tǒng)算法處理的效率[4]。本文基于 CUDA架構(gòu)對(duì)MTFC算法進(jìn)行并行化設(shè)計(jì)，并根據(jù)算法特征進(jìn)行優(yōu)化，提高了運(yùn)算速度，減少了算法執(zhí)行時(shí)間，滿足地面處理系統(tǒng)準(zhǔn)實(shí)時(shí)處理的要求。

2 M TFC技術(shù)現(xiàn)狀

基于MTFC的圖像復(fù)原是圖像增強(qiáng)處理技術(shù)中的一種，圖像復(fù)原前的遙感圖像是經(jīng)過成像鏈路退化和其間引入噪聲污染的圖像，圖像復(fù)原的目的就是處理這種圖像得到未被退化和污染的圖像的最佳估計(jì)。在衛(wèi)星地面處理系統(tǒng)中一般采用刀刃法進(jìn)行MTF在軌測量，并基于MTF曲線進(jìn)行復(fù)原，復(fù)原算法采用維納濾波[5-7]，MTFC過程如圖1所示，可以分為MTF曲線計(jì)算、MTF矩陣計(jì)算和圖像復(fù)原3個(gè)主要步驟，圖1中FFT和IFFT分別為傅里葉變換與傅里葉逆變換。

圖1 地面處理系統(tǒng)中的MTFC流程Fig.1 The flow of MTFC

在不考慮噪聲的情況下，刀刃法利用靶標(biāo)或圖像中低反射率目標(biāo)和高反射率目標(biāo)相交直線為參考，在直線實(shí)際亮度突變函數(shù)用f(x,y)表示，若相機(jī)的綜合響應(yīng)函數(shù)為h(x,y)，對(duì)于一個(gè)線性不變系統(tǒng)，則輸出圖像函數(shù)g(x,y)為[8]

式中 * 為卷積運(yùn)算符號(hào)；x,y為空間域像素坐標(biāo)；f(x,y)為原始圖像函數(shù)；g(x,y)為經(jīng)過大氣衰減、星載光學(xué)系統(tǒng)、傳感器、電子線路后退化的圖像函數(shù)；h(x,y)為綜合所有退化因素的函數(shù)，即點(diǎn)擴(kuò)散函數(shù)(Point Spread Function, PSF)。測量過程中，首先利用灰度進(jìn)行邊緣檢測，檢測出的邊緣點(diǎn)并非嚴(yán)格位于一條直線上，需要對(duì)邊界位置擬合、重采樣，然后進(jìn)行行配準(zhǔn)得到邊界位置的灰度分布，利用多項(xiàng)式擬合獲得邊界擴(kuò)散函數(shù)(Edge Spread Function, ESF)，對(duì)ESF進(jìn)行微分可以得到線擴(kuò)散函數(shù)(Line Spread Function,LSF)。對(duì)LSF進(jìn)行離散傅里葉變換，取變換之后各分量的模為各頻率的MTF值，并以第1個(gè)MTF值為基準(zhǔn)，作歸一化處理，就得到了一系列MTF值。由于截止頻率處的MTF值趨近于0，將頻率點(diǎn)以截至頻率為基準(zhǔn)作歸一化處理，則截止頻率為1，取頻率0～1處的MTF值構(gòu)成光學(xué)系統(tǒng)的MTF曲線。

建立MTFC模型必須要構(gòu)建二維的MTF矩陣。常規(guī)的處理方法是將水平MTF列向量乘以垂直MTF列向量，即：

式中 M TFu為在頻率u處水平的MTF值；M TFv為在頻率v處垂直的MTF值；M TF(u,v)為二維頻率坐標(biāo)為(u, v)處的MTF值。這種方法求得的45o方向的MTF值與水平或垂直方向的MTF值差別較大，為了消除這種差別，同時(shí)取水平與垂直方向0.5頻率處MTF值的平均值再衰減90%作為45o方向0.5頻率處的MTF值，再根據(jù)水平與垂直方向的MTF向量之間的比例關(guān)系進(jìn)行插值，即可得到二維插值MTF矩陣。由于模的對(duì)稱性，只需要求出 0～0.5頻率處的 MTF值，根據(jù)對(duì)稱性即可得到–0.5～0頻率處的MTF值，0.5頻率即截止頻率的一半。

基于 MTFC的圖像復(fù)原一般采用維納(Wiener)濾波算法，維納濾波復(fù)原是一種對(duì)噪聲起抑制和減少作用的方法[9]。維納濾波復(fù)原算法為

式中 F(u,v)為復(fù)原圖像的頻譜；G(u,v)為原始退化圖像的頻譜；Pf(u,v)和 Pn(u,v)分別為信號(hào)和噪聲的功率譜密度；k為一個(gè)可調(diào)系統(tǒng)，一般為0.02；* 表示共軛，即

式中 H(u,v)為二維MTF矩陣。維納濾波作為一種典型的有噪聲約束下輸入與估計(jì)最小均方誤差(MSE)的復(fù)原方法，可以有效解決在高頻處對(duì)噪聲的放大。

3 基于CUDA的MTFC算法

3.1 CUDA技術(shù)

CUDA是一種將GPU作為數(shù)據(jù)并行計(jì)算設(shè)備的軟硬件體系，其特點(diǎn)是將CPU作為主機(jī)，GPU作為協(xié)處理器或設(shè)備端，在這個(gè)模型中CPU和GPU協(xié)同工作，CPU負(fù)責(zé)邏輯性強(qiáng)的事務(wù)和串行計(jì)算任務(wù)，而GPU主要處理高度線程化的并行計(jì)算任務(wù)。在CUDA編程模型[10]中，如圖2所示，運(yùn)行在GPU上的CUDA并行計(jì)算函數(shù)稱為Kernel，一個(gè)完整的CUDA程序是由一系列的設(shè)備端Kernel函數(shù)并行步驟和主機(jī)端的串行處理步驟共同組成，而Kernel是整個(gè)CUDA程序中的一個(gè)可并行執(zhí)行步驟。在GPU端Kernel函數(shù)是以線程網(wǎng)絡(luò)(Grid)的形式組織，每個(gè)Grid再由若干個(gè)線程塊(Block)組成，每個(gè)Block中再包含很多個(gè)線程(Thread)。GPU線程的發(fā)起是輕量的, 其創(chuàng)建線程的系統(tǒng)開銷非常小, 線程切換所耗費(fèi)的時(shí)間也相當(dāng)短。

圖2 CUDA編程模型Fig.2 CUDA programm ing model

3.2 MTFC方法并行設(shè)計(jì)

基于MTFC的復(fù)原要求每次處理的圖像與二維MTF矩陣大小一致，而高分辨率、寬覆蓋影像比較大，需要對(duì)圖像進(jìn)行分塊處理。在地面系統(tǒng)處理中，為了消除維納濾波的邊緣影響，影像的四周邊緣不進(jìn)行復(fù)原處理，如圖3(a)所示，對(duì)abcd為頂點(diǎn)的圖像進(jìn)行處理時(shí)，上下左右各8個(gè)像元的邊緣保留原圖像，所以直接對(duì)圖像進(jìn)行分塊處理，會(huì)造成分塊邊界部分的復(fù)原質(zhì)量下降，為了保證復(fù)原質(zhì)量，本文設(shè)計(jì)了重疊分塊的方法，如圖3(b)所示，一個(gè)分塊與左及上一個(gè)分塊在行方向及列方向上都重疊16個(gè)像元，這樣就能保證復(fù)原算法在分塊邊緣位置的連續(xù)，從而得到與傳統(tǒng)算法一致的復(fù)原效果。

圖3 圖像并行分塊方法Fig.3 Image splitting method for parallel computiong

按照線程塊間并行和塊內(nèi)線程并行的結(jié)構(gòu)，本文設(shè)計(jì)了MTFC復(fù)原的算法，圖像復(fù)原的GPU算法流程如圖4所示。GPU程序一般是通過Kernel函數(shù)來實(shí)現(xiàn)的，本文定義5個(gè)主要順序執(zhí)行的Kernel函數(shù)：1) Kernel 1 將圖像量化類型轉(zhuǎn)換為適合傅里葉變換的復(fù)合類型；2) Kernel 2進(jìn)行傅里葉變換；3) Kernel 3 圖像濾波復(fù)原；4) Kernel 4 進(jìn)行傅里葉逆變換；5) Kernel 5 將復(fù)合類型轉(zhuǎn)為圖像量化類型。步驟中2)和4)采用CUDA中的CUFFT庫來實(shí)現(xiàn)。

3.3 并行優(yōu)化

直接將CPU程序在CUDA編程模型中進(jìn)行重構(gòu)，并不能帶來性能的顯著提升，只有對(duì)算法流程進(jìn)行充分分析，根據(jù)算法特征進(jìn)行設(shè)計(jì)和優(yōu)化，才能高效地利用GPU的計(jì)算資源，獲得數(shù)十倍乃至數(shù)百倍的加速比?；贑UDA實(shí)現(xiàn)并行算法，在數(shù)據(jù)分塊、內(nèi)核函數(shù)設(shè)計(jì)、存儲(chǔ)器使用等方面都有較多的選擇空間，因此基于CUDA實(shí)現(xiàn)并行算法可以有若干種實(shí)現(xiàn)方式，這樣算法的加速比、高效存儲(chǔ)器使用率、算法可以構(gòu)成一個(gè)三維空間。性能優(yōu)化的過程就是在三維空間逐漸尋找該最優(yōu)點(diǎn)的過程?；?CUDA實(shí)現(xiàn)的算法性能優(yōu)化主要從以下三方面考慮：1)最大限度的GPU端并行執(zhí)行；2)優(yōu)化存儲(chǔ)器使用以達(dá)到最大存儲(chǔ)器訪問帶寬；3)優(yōu)化指令流以獲取最大指令吞吐量[10]。

共享存儲(chǔ)器(Shared Memory)的合理使用可以提高線程的訪存速度。共享存儲(chǔ)器是GPU片內(nèi)的高速存儲(chǔ)器，它是可以被同一線程塊中所有線程訪問的可讀寫存儲(chǔ)器[11]，一般而言可以在一到兩個(gè)時(shí)鐘內(nèi)讀寫，因此使用共享存儲(chǔ)器取代全局存儲(chǔ)器會(huì)極大的節(jié)約帶寬。本文通過分析MTFC算法，將二維MTF矩陣轉(zhuǎn)化為一維矩陣，另外定義一個(gè)線程塊一次處理的圖像范圍為一個(gè)邏輯塊，將與該邏輯塊對(duì)應(yīng)大小的MTF值從GPU緩存拷入線程塊的共享存儲(chǔ)器，使得塊內(nèi)各線程通過共享的內(nèi)存來實(shí)現(xiàn)完全循環(huán)解纏處理，提高了處理速度。

在CUDA編程模型中，必須要用足夠多的活躍線程才能提高計(jì)算效能，隱藏訪存延時(shí)[12]。因此線程塊和邏輯塊的大小，會(huì)影響 MTFC算法的執(zhí)行效率，本文通過 CUDA效率分析器(CUDA Occupancy Calculator，COC)在GPU計(jì)算能力2.0的規(guī)格上進(jìn)行分析發(fā)現(xiàn)，當(dāng)線程塊大小為256(即16×16)，且每個(gè)線程塊使用的共享存儲(chǔ)器數(shù)據(jù)量為8 192字節(jié)(即邏輯塊大小為32×32)時(shí)，GPU線程束占用率達(dá)到最大值，如圖5、6所示，這種情況下，GPU的計(jì)算資源和存儲(chǔ)資源可以被完全使用。

圖4 圖像復(fù)原的GPU流程Fig.4 The flow of image restoration in GPU device

圖5 線程塊大小對(duì)效率的影響Fig.5 Impact of varying thread block size

圖6 共享存儲(chǔ)器大小對(duì)效率的影響Fig.6 Impact of varying shared memory usage

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

為驗(yàn)證算法的有效性，本文設(shè)計(jì)、編寫了基于CUDA的MTFC程序，并在硬件處理架構(gòu)上進(jìn)行實(shí)驗(yàn)。硬件平臺(tái)為Intel Core 2 Duo i7處理器，8G內(nèi)存，GPU平臺(tái)為NVIDIA Quadro 2000M。操作系統(tǒng)為Windows 7，編譯環(huán)境為M icrosoft Viusal Studio 2010、NVIDIA CUDA編譯器NVCC。實(shí)驗(yàn)數(shù)據(jù)為資源三號(hào)衛(wèi)星下視影像，圖像大小為24 530×24 575像素，16 bit量化，數(shù)據(jù)量為1.2 GB。圖7 為實(shí)驗(yàn)數(shù)據(jù)在CPU串行和GPU并行架構(gòu)下進(jìn)行MTFC處理的結(jié)果，由圖中可以看出，經(jīng)過復(fù)原處理，影像的邊緣比原圖更為清晰，細(xì)節(jié)更為豐富，通過圖像軟件比較，CPU串行MTFC處理及CUDA并行MTFC處理后的圖像結(jié)果完全一致。

圖7 待復(fù)原圖像和復(fù)原結(jié)果Fig.7 Image restoration results base on CPU and CUDA

為驗(yàn)證算法在不同數(shù)據(jù)量情況下的計(jì)算效率上，本文分別從實(shí)驗(yàn)數(shù)據(jù)截取 6 000×6 000、12 000×12 000、18 000×18 000、24 000×24 000像元作為形成模擬圖像，并對(duì)每個(gè)圖像運(yùn)行基于CUDA的MTFC程序和原CPU串行算法，連續(xù)運(yùn)行5次記錄運(yùn)行時(shí)間(不計(jì)文件IO時(shí)間)的平均值，表1列出了基于CUDA的MTFC算法與原基于CPU的MTFC算法的執(zhí)行時(shí)間和加速比。由表1中可以看出，基于CUDA的MTFC算法由于利用了GPU的眾核計(jì)算資源，使得算法的執(zhí)行時(shí)間大大縮短，圖像復(fù)原步驟可以獲得達(dá)到45倍的加速比，隨著圖像大小的增加，算法加速比逐漸增高，說明對(duì)于數(shù)據(jù)量越大的圖像，GPU加速效果越明顯，鑒于測試環(huán)境中GPU平臺(tái)只是NVIDIA公司用于圖形渲染的顯卡，如果移植到用于科學(xué)計(jì)算的Tesla系列專業(yè)GPU卡上，將會(huì)有更大程度的速度提升。

表1 不同圖像大小下的算法執(zhí)行平均時(shí)間對(duì)比Tab.1 Comparison of average executing time

另外，在選擇不同的線程塊和邏輯塊的大小時(shí)，MTFC算法的執(zhí)行時(shí)間會(huì)有差別，圖8列出了選擇不同的線程塊和邏輯塊大小進(jìn)行整景圖像MTFC算法執(zhí)行的時(shí)間。由圖8可以看出，在線程塊大小為16×16，邏輯塊大小為32×32的情況下，算法執(zhí)行時(shí)間最短，這與本文利用COC進(jìn)行并行優(yōu)化分析的結(jié)論是一致的。

圖8 不同線程塊和邏輯塊大小執(zhí)行時(shí)間對(duì)比Fig.8 Varing logical block size and thread block size

5 結(jié)束語

本文提出基于CUDA架構(gòu)的并行MTFC算法，通過分析MTFC算法流程及特征，采用線程塊間及塊內(nèi)線程間兩層并行模型，提高算法的并行度，另外，根據(jù)維納濾波復(fù)原算法的特點(diǎn)對(duì)存儲(chǔ)層次進(jìn)行優(yōu)化，提高共享存儲(chǔ)器的使用率。實(shí)驗(yàn)表明，該算法與原CPU串行算法處理結(jié)果質(zhì)量一致，但加速顯著，在NVIDIA一般顯卡環(huán)境下，即可達(dá)到45倍的加速比，大幅度提高了MTFC處理效率，如果移植到用于科學(xué)計(jì)算的Tesla系列專業(yè)GPU卡上，將會(huì)有更大程度的速度提升，在高分辨率衛(wèi)星地面處理系統(tǒng)中具有很好的應(yīng)用價(jià)值。

References)

[1]楊貴軍, 邢著榮. 小衛(wèi)星CCD相機(jī)MTF在軌測量與圖像復(fù)原[J].中國礦業(yè)大學(xué)學(xué)報(bào), 2011, 5(3): 481-486.YANG Guijun, XING Zhurong.On-orbit MTF Estimation and Restoration for CCD Cameras of Environment and Disaster Reduction Small Satellites(HJ1A/1B)[J]. Journal of China University of M ining & Technology, 2011, 5(3): 481-486. (in Chinese)

[2]Elsen E, Houstron M, Vishal V, et al. BN-body Simulation on GPUs[C]//Processing of 2006 ACM/IEEE Conference on Supercomputing, USA: New York, 2006: 188.

[3]Nvidia. CUDA Programm ing Guide Version 2.3.1[EB/OL]. http://developer.download.nvidia.com/compute/cuda/2_3/toolkit/docs/NVIDIA_CUDA_Programming_Guide_2.3.pdf, 2010-03-05.

[4]何國經(jīng), 劉德連.CUDA架構(gòu)下高光譜圖像光譜匹配的快速實(shí)現(xiàn)[J]. 航空兵器, 2011, 8(4): 3-6,12.HE Guojing,LIU Delian. High Speed Spectral Matching Approach for Hyperspectral Image Based on CUDA[J]. AERO WEAPONRY, 2011, 8(4): 3-6,12. (in Chinese)

[5]王先華, 喬延利, 易維寧.衛(wèi)星光學(xué)相機(jī)MTF在軌檢測方法研究[J]. 遙感學(xué)報(bào), 2007,11(3): 318-322.WANG Xianhua, QIAO Yanli, YI Weining. In-flight Measurement of Satellite Optical Camera MTF[J]. Journal of Remote Sensing, 2007,11(3): 318-322. (in Chinese)

[6]顧行發(fā), 李小英, 閔祥軍, 等. CBERS-02衛(wèi)星CCD相機(jī)MTF在軌測量及圖像MTF補(bǔ)償[J]. 中國科學(xué), 2005, 35(1): 26-39.GU Xingfa, LI Xiaoying, M IN Xiangjun, et al. In-flight MTF Measurement and MTFC of CBERS-02 Satellite CCD Images[J]. Science in China, 2005, 35(1): 26-39. (in Chinese)

[7]吳昀昭, 宮鵬. 北京一號(hào)小衛(wèi)星MTF 在軌測量與圖像復(fù)原[J]. 遙感應(yīng)用, 2007, (3): 49-53.WU Yunzhao, GONG Peng. On-orbit MTF Estimation and Restoration of Beijing-1 Images[J].Remote Sensing Application,2007, (3): 49-53. (in Chinese)

[8]周春平, 宮輝力, 李小娟, 等. 遙感圖像MTF復(fù)原國內(nèi)研究現(xiàn)狀[J]. 航天返回與遙感, 2009, 30(1): 27-32.ZHOU Chunping, GONG Huili, LI Xiaojuan, et al. The Summary of MTF Restoration on Remote Sensing Image[J].Spacecraft Recovery & Remote Sensing, 2009, 30(1): 27-32. (in Chinese)

[9]Kundur D, Hatzinakos D. Blind Image Deconvolution[J]. IEEE Signal Processing Magazine, 1996, 13(3): 43-64.

[10]胡婭. 基于GPU的BLAST程序的并行計(jì)算的研究[D]. 杭州: 浙江理工大學(xué), 2011.HU Ya. Research on GPU-based Parallel Computing on BLAST Program[D]. Hangzhou: Zhejiang Sci-Tech University, 2011.(in Chinese)

[11]張舒, 褚艷利. GPU高性能運(yùn)算之CUDA[M]. 北京: 中國水利水電出版社, 2009: 46.ZHANG Shu, ZHU Yanli. GPU High Performance Computing –CUDA[M]. Beijing:China Water Power Press, 2009: 46. (in Chinese)

[12]Muthu M B, Rajesh B. Optim izing Sparse Matrix-vector Multipli-cation on GPUs[R]. [S.1.]: IBM, 2008.