方留楊，王 密，李德仁，潘 俊

武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079

負(fù)載分配的CPU／GPU高分辨率衛(wèi)星影像調(diào)制傳遞補(bǔ)償方法

方留楊，王 密，李德仁，潘 俊

武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079

系統(tǒng)地探討了使用CPU／GPU協(xié)同處理理論對(duì)高分辨率衛(wèi)星影像進(jìn)行MTF補(bǔ)償?shù)姆椒?。首先在GPU上對(duì)方法進(jìn)行了基本實(shí)現(xiàn)，并通過(guò)三種性能優(yōu)化策略（執(zhí)行配置優(yōu)化、存儲(chǔ)訪問(wèn)優(yōu)化和指令優(yōu)化）進(jìn)一步提高了方法的執(zhí)行效率。在Intel Xeon E5650 CPU和NVIDIA Tesla C2050 GPU組成的CPU／GPU系統(tǒng)中對(duì)高分一號(hào)衛(wèi)星全色影像進(jìn)行MTF補(bǔ)償，加速比達(dá)到42．80倍。在此基礎(chǔ)上，為充分利用CPU的計(jì)算性能，使用CPU／GPU負(fù)載分配策略將部分負(fù)載分配給CPU進(jìn)行處理，使用該策略后，方法加速比達(dá)到47．82倍，相應(yīng)的處理時(shí)間壓縮至1．62 s，可滿(mǎn)足對(duì)高分辨率衛(wèi)星影像進(jìn)行近實(shí)時(shí)MTF補(bǔ)償?shù)男枨蟆?/p>

CPU／GPU；調(diào)制傳遞函數(shù)補(bǔ)償；執(zhí)行配置優(yōu)化；存儲(chǔ)訪問(wèn)優(yōu)化；指令優(yōu)化；負(fù)載分配策略

1 引 言

高分辨率衛(wèi)星影像在成像過(guò)程中受大氣擾動(dòng)、平臺(tái)震顫、傳感器性能下降等因素影響，致使影像質(zhì)量退化。基于調(diào)制傳遞函數(shù)（modulation transfer function，MTF）理論的圖像復(fù)原方法是目前較為常用的改善高分辨率衛(wèi)星影像質(zhì)量的方法。該方法通過(guò)計(jì)算成像系統(tǒng)的MTF曲線在不同空間頻段的下降程度來(lái)反推實(shí)際影像頻譜中高頻部分的上升程度，達(dá)到改善影像質(zhì)量的目的［1］。由于需要對(duì)影像進(jìn)行較為復(fù)雜的逐像素頻域?yàn)V波處理，因此MTF補(bǔ)償?shù)挠?jì)算量很大。特別是在衛(wèi)星數(shù)據(jù)源快速增多、空間分辨率大幅提高、需要處理的影像數(shù)據(jù)量急劇增長(zhǎng)的背景下，傳統(tǒng)的MTF補(bǔ)償計(jì)算時(shí)間長(zhǎng)、效率低下等問(wèn)題日益明顯，已無(wú)法滿(mǎn)足在許多應(yīng)急應(yīng)用領(lǐng)域（如軍事情報(bào)判讀、搶險(xiǎn)救災(zāi)）中對(duì)衛(wèi)星影像進(jìn)行實(shí)時(shí)或準(zhǔn)實(shí)時(shí)MTF補(bǔ)償處理的要求，成為制約整個(gè)高分辨率衛(wèi)星數(shù)據(jù)處理快速完成的瓶頸。

近年來(lái)，隨著半導(dǎo)體技術(shù)和計(jì)算機(jī)技術(shù)的發(fā)展，圖形處理單元（graphics processing unit，GPU）的計(jì)算性能以超過(guò)摩爾定律的速度不斷提高，已成為當(dāng)前計(jì)算機(jī)系統(tǒng)中具備高性能處理能力的部件。目前，國(guó)內(nèi)外相關(guān)學(xué)者對(duì)使用GPU提高計(jì)算密集型遙感數(shù)據(jù)處理算法的效率開(kāi)展了一些研究：文獻(xiàn)［2］使用GPU進(jìn)行高光譜影像非監(jiān)督波段選取試驗(yàn)，算法效率提高22倍；文獻(xiàn)［3］使用GPU對(duì)光學(xué)航空影像進(jìn)行幾何校正，算法加速比達(dá)到45倍；文獻(xiàn)［4］在GPU上建立紅外大氣探測(cè)干涉儀的輻射傳輸模型，獲得了上百倍的加速比。在國(guó)內(nèi)，文獻(xiàn)［5］實(shí)現(xiàn)了基于GPU的并行影像匹配算法，將影像匹配速度提高了7倍。從上述應(yīng)用可以看出，雖然各算法加速比不一致（與硬件環(huán)境、算法特點(diǎn)有關(guān)），但處理效率都明顯提高。MTF補(bǔ)償計(jì)算密集，且逐像元進(jìn)行，具有先天并行性，非常適合于在GPU上進(jìn)行處理。文獻(xiàn)［6］使用CPU／GPU對(duì)資源三號(hào)衛(wèi)星影像進(jìn)行MTF補(bǔ)償，取得了較好的效果，但其性能優(yōu)化并不充分，算法效率仍有進(jìn)一步提高空間。本文使用基于Fermi架構(gòu)的新一代GPU，綜合考慮充分利用CPU的計(jì)算性能，系統(tǒng)地探討使用CPU／GPU協(xié)同處理理論對(duì)高分辨率衛(wèi)星影像進(jìn)行MTF補(bǔ)償?shù)姆椒?，涉及方法在GPU上的實(shí)現(xiàn)、性能的充分優(yōu)化以及CPU和GPU間負(fù)載分配等一系列問(wèn)題。

2 原理與方法

2．1 MTF補(bǔ)償原理

設(shè)f（x，y）為原始未退化影像，g（x，y）為觀測(cè)到的退化影像，退化的數(shù)學(xué)模型為［7］

式中，h（x，y）為點(diǎn)擴(kuò)散函數(shù)；＊代表卷積操作。

對(duì)兩邊進(jìn)行傅里葉變換，得頻域表達(dá)形式

式中，MTF（u，v）是二維MTF矩陣；kw是與圖像信噪比有關(guān)的先驗(yàn)常數(shù)。

2．2 CPU／GPU協(xié)同處理

在當(dāng)前主流的CPU／GPU協(xié)同處理系統(tǒng)中，通常包含一個(gè)多核CPU和一片眾核GPU。其中在多核CPU 上應(yīng)用最廣泛的編程模型是Open MP。Open MP是共享存儲(chǔ)系統(tǒng)并行編程的工業(yè)標(biāo)準(zhǔn)，提供對(duì)并行區(qū)域、工作區(qū)間共享、同步等處理的支持。Open MP可直接在高級(jí)語(yǔ)言（如C、Fortran）的串行代碼上實(shí)現(xiàn)并行：在源代碼中嵌入編譯制導(dǎo)語(yǔ)句，編譯后執(zhí)行時(shí)制導(dǎo)語(yǔ)句啟動(dòng)多線程環(huán)境實(shí)現(xiàn)并行計(jì)算，不需要對(duì)程序做大的修改［13］。關(guān)于Open MP的詳細(xì)內(nèi)容可參見(jiàn)文獻(xiàn)［13］。

早期的GPU開(kāi)發(fā)直接使用圖形學(xué)API（如Open GL，Direct X等）進(jìn)行編程，程序員需要對(duì)圖形學(xué)硬件和編程接口有深入的了解，開(kāi)發(fā)難度較大，因此應(yīng)用并不廣泛。2007年，NVIDIA公司推出了統(tǒng)一計(jì)算設(shè)備架構(gòu)（compute unified device architecture，CUDA），該架構(gòu)通過(guò)簡(jiǎn)單易式中，H（u，v）的值是復(fù)數(shù)，若對(duì)幅值作歸一化，使得零頻率的幅值為1，則稱(chēng)此歸一化的幅值為調(diào)制傳遞函數(shù)（MTF）［1］。因此，退化模型可變?yōu)?/p>

計(jì)算上式中MTF的理論和方法有許多，主要包括點(diǎn)脈沖法、正弦輸入法、刀刃法和脈沖法［810］，其中后兩種方法運(yùn)用情況更為普遍。這兩種方法需要從衛(wèi)星在軌運(yùn)行獲取的影像中摳取靶標(biāo)子圖，通過(guò)子圖中的邊緣地物或線性地物生成MTF曲線，并構(gòu)造二維MTF矩陣。由于靶標(biāo)子圖較小，計(jì)算量不大，因此所需時(shí)間較短。關(guān)于MTF曲線生成和二維MTF矩陣構(gòu)造可參見(jiàn)文獻(xiàn)［10—11］。

MTF補(bǔ)償是退化模型的逆過(guò)程，在實(shí)際應(yīng)用中對(duì)退化影像進(jìn)行MTF補(bǔ)償時(shí)，需要選取一個(gè)包含MTF矩陣的濾波算子P，其表達(dá)形式為

常用的濾波算子包括逆濾波、維納濾波、修正反轉(zhuǎn)濾波、約束最小二乘方濾波和Lucy－Richardson濾波等。其中維納濾波考慮了噪聲的影響，是一種很常用的圖像恢復(fù)算子，可直接應(yīng)用于頻率域［7，1112］。頻率域的維納濾波算子如下式所示用的類(lèi)C語(yǔ)言進(jìn)行開(kāi)發(fā)，能夠比較容易地獲得幾倍、幾十倍，乃至上百倍的加速效率，得到了非常廣泛的使用［14］。CUDA實(shí)際上是一種“串行—并行”混合的CPU／GPU協(xié)同編程模型：需要并行處理的部分被組織成“核函數(shù)”，“核函數(shù)”啟動(dòng)時(shí)產(chǎn)生大量線程，這些線程以線程網(wǎng)格和線程塊的形式進(jìn)行層次性組織，同時(shí)對(duì)數(shù)據(jù)進(jìn)行并行化處理；其余串行執(zhí)行步驟由CPU進(jìn)行處理，這些步驟主要包括在CPU和GPU間傳輸數(shù)據(jù)、配置核函數(shù)執(zhí)行參數(shù)、啟動(dòng)核函數(shù)并行執(zhí)行程序等［1517］。此外，CUDA設(shè)計(jì)了一個(gè)層次性存儲(chǔ)體系，包括全局存儲(chǔ)器、常數(shù)存儲(chǔ)器、共享存儲(chǔ)器等，這些存儲(chǔ)器的存儲(chǔ)容量和訪問(wèn)速度互不相同。關(guān)于CUDA編程模型、線程組織結(jié)構(gòu)和層次性存儲(chǔ)體系等詳細(xì)內(nèi)容可參見(jiàn)文獻(xiàn)［15，18］。

2．3 負(fù)載分配的CPU／GPU MTF補(bǔ)償方法

在CPU／GPU協(xié)同處理系統(tǒng)中進(jìn)行MTF補(bǔ)償時(shí)，需要單獨(dú)指派一個(gè)CPU核與GPU進(jìn)行交互（CPU和GPU之間的數(shù)據(jù)傳輸、核函數(shù)執(zhí)行參數(shù)配置、核函數(shù)啟動(dòng)等），本文將其命名為“控制核”；此外，為充分利用系統(tǒng)中其余CPU核的計(jì)算性能，將其作為“計(jì)算核”，與GPU一起共同完成對(duì)負(fù)載（即退化影像）的MTF補(bǔ)償處理，如圖1所示。其中GPU方法實(shí)現(xiàn)及性能優(yōu)化、控制核和GPU之間的交互使用CUDA實(shí)現(xiàn)；計(jì)算核并行區(qū)域構(gòu)造、控制核與計(jì)算核之間的交互使用Open MP實(shí)現(xiàn)。GPU和計(jì)算核之間的負(fù)載分配依據(jù)負(fù)載分配策略進(jìn)行，該策略根據(jù)計(jì)算核的數(shù)量和方法的實(shí)際加速比為計(jì)算核和GPU合理分配對(duì)應(yīng)的負(fù)載量，使得方法的總執(zhí)行效率達(dá)到最優(yōu)。

圖1 負(fù)載分配的CPU／GPU MTF補(bǔ)償方法Fig．1 The workload－distribution based CPU／GPU MTF compensation approach

2．3．1 GPU方法實(shí)現(xiàn)及性能優(yōu)化

由于需要進(jìn)行逐像素頻域?yàn)V波，計(jì)算較為復(fù)雜，因此GPU上的MTF補(bǔ)償被劃分為3個(gè)核函數(shù)實(shí)現(xiàn)：Convert ToComplex，Compensation和ConvertFromComplex。其中Convert ToComplex將影像數(shù)據(jù)類(lèi)型轉(zhuǎn)為在GPU上進(jìn)行傅里葉變換所需的Complex數(shù)據(jù)類(lèi)型，Compensation執(zhí)行頻域維納濾波，Convert From Complex將Complex數(shù)據(jù)類(lèi)型轉(zhuǎn)回為影像數(shù)據(jù)類(lèi)型。傅里葉變換使用NVIDIA提供的CUDA快速傅里葉變換（CUDA fast Fourier transform，CUFFT）函數(shù)庫(kù)實(shí)現(xiàn)。在核函數(shù)中，將線程塊中線程數(shù)量設(shè)置為Fermi架構(gòu)GPU的上限1024（32×32），且每個(gè)線程負(fù)責(zé)對(duì)一個(gè)像素進(jìn)行處理。因此，若設(shè)退化影像的大小為lwidth×lheight，則共需產(chǎn)生［lwidth／32］×［lheight／32］個(gè)線程塊完成對(duì)整幅影像的處理（進(jìn)行向上取整的原因是保證所有像素皆被覆蓋）。

上述方法僅為基本實(shí)現(xiàn)，可保證MTF補(bǔ)償在GPU上順利完成，但無(wú)法充分發(fā)揮GPU的計(jì)算能力。因此需對(duì)方法進(jìn)行性能優(yōu)化，主要通過(guò)以下3個(gè)方面展開(kāi)：

（1）執(zhí)行配置優(yōu)化可進(jìn)一步細(xì)分為提高流式多處理器中線程占有率和線程多元素重訪兩個(gè)策略。在提高流式多處理器中線程占有率方面，基于Fermi架構(gòu)的GPU需要在滿(mǎn)足如下兩個(gè)限制條件的前提下，使流式多處理器中線程占有率達(dá)到最高（流式多處理器中包含1536個(gè)線程）［18］，以提高GPU線程調(diào)度的性能：①線程塊中線程數(shù)量小于1024；②流式多處理器中線程塊數(shù)量小于8。

經(jīng)過(guò)分析，當(dāng)線程塊中線程數(shù)量設(shè)置為256或512時(shí)，流式多處理器中線程塊數(shù)量分別為6或3，同時(shí)滿(mǎn)足限制條件。在這兩種情況下，流式多處理器都包含1536個(gè)線程，線程占有率達(dá)到最高。

線程多元素重訪策略利用GPU上多次顯存訪問(wèn)可被流水化以及線程尋址可重用等特點(diǎn)，令一個(gè)GPU線程同時(shí)處理多個(gè)像素，以提高GPU上線程運(yùn)行的性能。需要注意的是，GPU線程運(yùn)行的性能并不會(huì)隨著線程訪問(wèn)像素量的增加而無(wú)限增長(zhǎng)，而是存在一個(gè)峰值，因此需根據(jù)實(shí)際情況找出性能達(dá)到峰值時(shí)GPU線程對(duì)應(yīng)處理的像素量，使算法執(zhí)行效率達(dá)到最優(yōu)。

（2）存儲(chǔ)訪問(wèn)優(yōu)化也可細(xì)分為兩個(gè)策略：存儲(chǔ)層次性訪問(wèn)和傳輸計(jì)算堆疊。在存儲(chǔ)層次性訪問(wèn)方面，基于Fermi架構(gòu)的GPU除全局存儲(chǔ)器外，還提供了由共享存儲(chǔ)器和一級(jí)緩存共同組成的容量為64 KB的片上緩存。該片上緩存可柔性重組，提供兩種容量分配模式：①48 KB共享存儲(chǔ)器＋16 KB一級(jí)緩存；②16 KB共享存儲(chǔ)器＋48 KB一級(jí)緩存［19］。在核函數(shù)Compensation進(jìn)行維納濾波的過(guò)程中，MTF矩陣被多次使用，因此可將其載入至共享存儲(chǔ)器中，減少程序訪問(wèn)全局存儲(chǔ)器的次數(shù)；又為了一次性盡可能多地載入MTF矩陣，故選擇分配模式①，并通過(guò)式（6）確定每個(gè)線程對(duì)應(yīng)處理的像素量，使得共享存儲(chǔ)器空間被完全利用。Convert ToComplex和Convert FromComplex兩個(gè)核函數(shù)無(wú)數(shù)據(jù)重復(fù)性訪問(wèn)，因此選擇分配模式②，將更多影像緩存至一級(jí)cache；由于沒(méi)有使用共享存儲(chǔ)器，因此這兩個(gè)核函數(shù)的每個(gè)線程仍然對(duì)應(yīng)處理一個(gè)像素。

每線程對(duì)應(yīng)處理像素量＝

在CUDA中，數(shù)據(jù)在內(nèi)存和GPU顯存之間的傳輸默認(rèn)為串行阻塞式傳輸，即數(shù)據(jù)完全從內(nèi)存拷貝至GPU顯存后，GPU核函數(shù)才能執(zhí)行；當(dāng)GPU核函數(shù)執(zhí)行完畢后，數(shù)據(jù)再?gòu)腉PU顯存中拷回至內(nèi)存。基于Fermi架構(gòu)的GPU內(nèi)部固化有一個(gè)數(shù)據(jù)拷貝引擎和一個(gè)核函數(shù)執(zhí)行引擎，因此可使用流機(jī)制來(lái)實(shí)現(xiàn)數(shù)據(jù)傳輸和計(jì)算的堆疊。流是一個(gè)“先進(jìn)先出”的命令流水線［18］，包含1次內(nèi)存至顯存數(shù)據(jù)拷貝操作、5次核函數(shù)執(zhí)行操作（Convert ToComplex、傅里葉正變換、Compensation、傅里葉逆變換和Convert From Complex）和1次顯存至內(nèi)存數(shù)據(jù)拷貝操作。當(dāng)程序開(kāi)始執(zhí)行時(shí)，GPU底層硬件將數(shù)據(jù)拷貝操作映射至數(shù)據(jù)拷貝引擎，將核函數(shù)執(zhí)行操作映射至核函數(shù)執(zhí)行引擎，隨后這兩個(gè)引擎異步執(zhí)行相應(yīng)的任務(wù)，實(shí)現(xiàn)數(shù)據(jù)傳輸和計(jì)算的堆疊。

（3）指令優(yōu)化。為達(dá)到更高的處理效率，GPU上放置了數(shù)百個(gè)計(jì)算核心，遠(yuǎn)多于CPU上計(jì)算核心的數(shù)量。因此與CPU相比，GPU上控制部件較少，使得GPU對(duì)控制邏輯的處理能力較弱：一些高級(jí)語(yǔ)言中常用的控制邏輯（如線程束內(nèi)部“if－else”分支、“for”循環(huán)等）在GPU上的執(zhí)行時(shí)間遠(yuǎn)長(zhǎng)于在CPU上的執(zhí)行時(shí)間，導(dǎo)致GPU性能損失。在3個(gè)核函數(shù)中，由于橫縱坐標(biāo)方向的線程數(shù)大于影像像素個(gè)數(shù)，因此在邊界像素處理時(shí)引入了邏輯分支，但這些分支語(yǔ)句為逐像素分支，不屬于線程束內(nèi)部的邏輯分支，故不會(huì)導(dǎo)致性能損失；在核函數(shù)Compensation中，為充分利用共享存儲(chǔ)器，每個(gè)線程對(duì)應(yīng)處理的像素量會(huì)多于1個(gè)，因此會(huì)產(chǎn)生“for”循環(huán)的使用，并使程序執(zhí)行性能下降；為提高性能，需進(jìn)行循環(huán)解纏處理（將循環(huán)邏輯展開(kāi)為順序邏輯執(zhí)行）。Convert－ToComplex和Convert FromComplex兩個(gè)核函數(shù)中每個(gè)線程仍對(duì)應(yīng)一個(gè)像素的工作量，因此無(wú)需循環(huán)解纏。循環(huán)解纏前后核函數(shù)Compensation的執(zhí)行流程如表1所示（假設(shè)每個(gè)線程對(duì)應(yīng)處理n個(gè)像素）。

表1 循環(huán)解纏前后核函數(shù)Compensation的執(zhí)行流程Tab．1 The original execution flow of kernel“Compensation”and its loop unroll counterpart

綜上所述，性能優(yōu)化后的CPU／GPU MTF補(bǔ)償流程如圖2所示，其中灰色標(biāo)記的步驟表示使用流實(shí)現(xiàn)了傳輸計(jì)算堆疊。

圖2 性能優(yōu)化后的CPU／GPU MTF補(bǔ)償流程（灰色標(biāo)記步驟表示使用流實(shí)現(xiàn)了傳輸計(jì)算堆疊）Fig．2 The CPU／GPU MTF compensation approach after performance tuning（the gray－marked steps are streamed for the overlap of data transfer and computation）

2．3．2 CPU／GPU負(fù)載分配策略

為充分利用系統(tǒng)中計(jì)算核的計(jì)算性能，將其與GPU一起共同完成對(duì)負(fù)載（即退化影像）的MTF補(bǔ)償處理。在這種情況下，需為兩者合理分配負(fù)載量，使方法的總執(zhí)行效率達(dá)到最優(yōu)。令P計(jì)、PGPU、m和α分別代表計(jì)算核的理論峰值計(jì)算性能（假設(shè)所有計(jì)算核的計(jì)算性能都相同）、GPU的理論峰值計(jì)算性能、計(jì)算核的數(shù)量以及分配給GPU的負(fù)載率。若總負(fù)載量為L(zhǎng)，則方法總執(zhí)行時(shí)間T由下式確定

容易發(fā)現(xiàn)，當(dāng)計(jì)算核和GPU同時(shí)完成處理時(shí)（T計(jì)＝TGPU），整個(gè)方法執(zhí)行時(shí)間最短，即當(dāng)α由下式?jīng)Q定時(shí)，T值最小

通常情況下，計(jì)算核和GPU的理論峰值計(jì)算性能都無(wú)法達(dá)到。因此，在實(shí)際應(yīng)用中，可用MTF補(bǔ)償?shù)腉PU版本相對(duì)CPU串行版本的加速比替代PGPU和P計(jì)。令s代表加速比，替換PGPU／P計(jì)，則由式（8）可得

3 試驗(yàn)結(jié)果與分析

2013年4月26日，我國(guó)高分辨率對(duì)地觀測(cè)衛(wèi)星系統(tǒng)首星“高分一號(hào)”成功發(fā)射升空，該衛(wèi)星搭載兩臺(tái)2 m分辨率全色／8 m分辨率多光譜相機(jī)和4臺(tái)16 m分辨率多光譜相機(jī)。其中全色相機(jī)采用10 bit量化，標(biāo)準(zhǔn)景影像大小為18 192像素×18 499像素，數(shù)據(jù)量較大（約642 MB），具有較好的代表性，因此本節(jié)選擇該影像作為試驗(yàn)影像。

本節(jié)使用加速比來(lái)評(píng)價(jià)本文方法相對(duì)于傳統(tǒng)串行方法執(zhí)行效率的提升，定義如下

式中，TCPU為傳統(tǒng)串行方法在CPU上的執(zhí)行時(shí)間；TCPU／GPU為本文方法在CPU／GPU系統(tǒng)上的執(zhí)行時(shí)間，由兩部分組成：數(shù)據(jù)在內(nèi)存和GPU顯存之間的傳輸時(shí)間和GPU實(shí)際執(zhí)行時(shí)間；s為加速比。

此外，使用性能提升比來(lái)評(píng)價(jià)某一優(yōu)化手段作用前后方法執(zhí)行效率的提升，定義如下

式中，Tbef和Taft分別為某一優(yōu)化手段作用前后方法的執(zhí)行時(shí)間；p為性能提升比。

3．1 方法性能

本文試驗(yàn)環(huán)境中的CPU為Intel Xeon E5650，含6個(gè)處理核心，單精度浮點(diǎn)峰值運(yùn)算能力約為128 Gflops；GPU選擇當(dāng)前主流的基于Fermi架構(gòu)的NVIDIA Tesla C2050，含448個(gè)處理核心，單精度浮點(diǎn)峰值運(yùn)算能力約為1030 Gflops。此外，NVIDIA Tesla C2050的片上緩存大小為64 KB（可柔性重組），全局存儲(chǔ)器容量為3 GB，內(nèi)部含有一個(gè)數(shù)據(jù)拷貝引擎和一個(gè)核函數(shù)執(zhí)行引擎；主機(jī)內(nèi)存大小為24 GB；全局存儲(chǔ)器和內(nèi)存之間通過(guò)PCIE2．0×16接口進(jìn)行通信；試驗(yàn)平臺(tái)為Win32（在此平臺(tái)上編譯和執(zhí)行的程序?yàn)?2位）。由于影像數(shù)據(jù)量較大，超過(guò)了32位程序一次能分配的最大內(nèi)存，本試驗(yàn)將影像劃分為1024像素×1024像素大小的數(shù)據(jù)塊進(jìn)行分塊載入和處理，保證內(nèi)存的成功分配。試驗(yàn)影像標(biāo)準(zhǔn)景大小為18 192像素×18 499像素，若將左上邊數(shù)據(jù)塊與影像左上角對(duì)齊，則劃分給最右邊和最下邊數(shù)據(jù)塊的像素不足1024（分別為784和67），可采取將無(wú)像素部分填0的方法將數(shù)據(jù)塊補(bǔ)齊至1024像素×1024像素，然后再與其余數(shù)據(jù)塊相同的方式進(jìn)行MTF補(bǔ)償，補(bǔ)償完畢后只保存對(duì)應(yīng)像素即可。

首先測(cè)得在試驗(yàn)環(huán)境的CPU上MTF曲線生成、二維MTF矩陣構(gòu)造以及MTF補(bǔ)償?shù)膱?zhí)行時(shí)間，如表2所示。其中MTF曲線生成選用刀刃法，作用于從退化影像扣取的靶標(biāo)子圖上（30像素×30像素，1．8 KB）；二維MTF矩陣構(gòu)造選擇文獻(xiàn)［11］中方法，矩陣大小與數(shù)據(jù)塊大小保持一致（1024像素×1024像素，4 MB）；MTF補(bǔ)償按分塊策略進(jìn)行，作用于高分一號(hào)全色標(biāo)準(zhǔn)景影像上（18 192像素×18 499像素，641 MB）。補(bǔ)償過(guò)程中的傅里葉正／反變換采用快速傅里葉變換（fast Fourier transform in the west，F(xiàn)FTW）函數(shù)庫(kù)實(shí)現(xiàn)。

表2 CPU上MTF曲線生成、二維MTF矩陣構(gòu)造和MTF補(bǔ)償執(zhí)行時(shí)間Tab．2 The run time of MTF curve generation，2D MTF matrix construction and MTF compensation on CPU s

從表中可知，與MTF曲線生成和二維MTF矩陣構(gòu)造相比，MTF補(bǔ)償執(zhí)行時(shí)間最長(zhǎng)，約占整個(gè)圖像復(fù)原處理執(zhí)行時(shí)間的99．7%，成為制約整個(gè)處理流程快速完成的瓶頸。故需使用本文方法對(duì)其進(jìn)行加速處理，提高處理效率。

經(jīng)測(cè)試，在試驗(yàn)環(huán)境的GPU上，僅進(jìn)行基本實(shí)現(xiàn)（未進(jìn)行性能優(yōu)化）的MTF補(bǔ)償方法的執(zhí)行時(shí)間為3．39 s，對(duì)應(yīng)的加速比達(dá)到22．85倍。接下來(lái)使用3種性能優(yōu)化策略對(duì)方法進(jìn)行優(yōu)化：首先進(jìn)行執(zhí)行配置優(yōu)化，將線程塊中線程數(shù)量設(shè)置為256，使流式多處理器中的線程占有率達(dá)到最高，并使用線程多元素重訪策略令一個(gè)GPU線程處理多個(gè)像素，圖3列出了隨著GPU線程處理像素?cái)?shù)量變化，方法性能提升比的變化情況，從圖中可以發(fā)現(xiàn)，當(dāng)GPU線程對(duì)應(yīng)處理4個(gè)像素時(shí)，方法的執(zhí)行效率最優(yōu)，對(duì)應(yīng)的性能提升比為5．1%；其次使用存儲(chǔ)層次性訪問(wèn)策略，對(duì)片上緩存進(jìn)行動(dòng)態(tài)分配并將MTF矩陣載入至共享存儲(chǔ)器中，又為了充分利用共享存儲(chǔ)器，根據(jù)式（6）計(jì)算核函數(shù)Compensation中每個(gè)線程對(duì)應(yīng)處理的像素量：動(dòng)態(tài)分配后共享存儲(chǔ)器容量為48 KB；流式多處理器中線程數(shù)量經(jīng)優(yōu)化后為1536；數(shù)據(jù)類(lèi)型為4字節(jié)單精度浮點(diǎn)型（Complex數(shù)據(jù)類(lèi)型中的實(shí)部［20］），因此可算出每個(gè)線程對(duì)應(yīng)8個(gè)像素的工作量；接下來(lái)使用流實(shí)現(xiàn)傳輸計(jì)算堆疊，從圖4可知，當(dāng)創(chuàng)建4個(gè)流的時(shí)候，方法的性能提升比最高（6．5%），由于流的創(chuàng)建和銷(xiāo)毀會(huì)帶來(lái)額外的開(kāi)銷(xiāo)，因此創(chuàng)建過(guò)多的流反而會(huì)導(dǎo)致方法性能下降；最后進(jìn)行指令優(yōu)化，對(duì)核函數(shù)Compensation進(jìn)行循環(huán)解纏處理。

圖3 使用線程多元素重訪策略后方法的性能提升比Fig．3 The performance improvement ratio after assigning each GPU thread the workload of multiple pixels

圖5列出了逐步使用3種性能優(yōu)化策略后，在試驗(yàn)環(huán)境的GPU上MTF補(bǔ)償?shù)膱?zhí)行時(shí)間和加速比。可以看出，執(zhí)行配置優(yōu)化后，方法加速比從22．85倍提高到33．11倍，性能提升最為明顯，達(dá)到44．9%；在此基礎(chǔ)上進(jìn)行存儲(chǔ)訪問(wèn)優(yōu)化后，方法性能進(jìn)一步提升，加速比達(dá)到42．10倍；當(dāng)對(duì)核函數(shù)Compensation進(jìn)行循環(huán)解纏處理后，整個(gè)方法的處理時(shí)間壓縮至1．81 s，相應(yīng)的加速比達(dá)到42．80倍。此外，還采用了文獻(xiàn)［6］中的方法對(duì)試驗(yàn)影像進(jìn)行了MTF補(bǔ)償，由于未使用線程多元素重訪和傳輸計(jì)算堆疊進(jìn)行性能優(yōu)化，整個(gè)方法的執(zhí)行時(shí)間為2．02 s，相應(yīng)的加速比僅為38．35倍。故本文方法的執(zhí)行效率較優(yōu)。

圖4 使用流實(shí)現(xiàn)傳輸計(jì)算堆疊后方法的性能提升比Fig．4 The performance improvement ratio after creating multiple streams

圖5 逐步使用性能優(yōu)化策略后GPU上MTF補(bǔ)償?shù)膱?zhí)行時(shí)間和加速比Fig．5 The run time and speedup ratio of MTF Compensation after gradually applying three optimization approaches on GPU

接下來(lái)使用CPU／GPU負(fù)載分配策略將部分負(fù)載（退化影像）分配給CPU計(jì)算核進(jìn)行處理，進(jìn)一步提高方法的執(zhí)行效率。試驗(yàn)環(huán)境中的CPU為6核Intel Xeon E5650，其中1個(gè)核用作控制核，其余5個(gè)核可作為計(jì)算核；此外，方法加速比已知。因此可通過(guò)式（9）計(jì)算出理論上分配給GPU和計(jì)算核的負(fù)載率。令式（9）中m為5，s為42．80，可得到α的理論值為0．895。即理論上，分配給GPU的負(fù)載百分比為89．5%；相應(yīng)地，分配給計(jì)算核的負(fù)載百分比為10．5%。在具體實(shí)施時(shí)，使用Open MP的parallel子句創(chuàng)建6個(gè)線程，然后通過(guò)相應(yīng)的Open MP API（omp＿get＿thread＿num（）函數(shù)）獲取線程號(hào)并給其指派任務(wù)，其中1個(gè)線程作為控制核線程與GPU交互，其余5個(gè)線程作為計(jì)算核線程，共同對(duì)分配給計(jì)算核的負(fù)載進(jìn)行并行處理（其中正反傅里葉變換的并行化通過(guò)調(diào)用FFTW函數(shù)庫(kù)內(nèi)嵌的多線程接口實(shí)現(xiàn)，影像數(shù)據(jù)類(lèi)型和Complex數(shù)據(jù)類(lèi)型的相互轉(zhuǎn)換以及維納濾波的并行化通過(guò)線程間任務(wù)分配進(jìn)行）。

為驗(yàn)證理論值的正確性，本文實(shí)測(cè)了在GPU對(duì)應(yīng)負(fù)載率變化的情況下GPU和CPU的執(zhí)行時(shí)間，如圖6所示。從圖中發(fā)現(xiàn)，當(dāng)GPU對(duì)應(yīng)負(fù)載率為0．893時(shí)，GPU和CPU的執(zhí)行時(shí)間相等（1．62 s），實(shí)際值與理論值之間存在一個(gè)大小為0．002細(xì)微誤差?？蓪⒃撜`差歸結(jié)于CUDA運(yùn)行時(shí)環(huán)境初始化帶來(lái)的時(shí)間開(kāi)銷(xiāo)，使得GPU執(zhí)行時(shí)間不與其對(duì)應(yīng)負(fù)載率嚴(yán)格成正比關(guān)系：當(dāng)GPU對(duì)應(yīng)負(fù)載率趨近于0時(shí)，GPU執(zhí)行時(shí)間并不趨近于0（當(dāng)GPU對(duì)應(yīng)負(fù)載率僅為0．005時(shí)，GPU執(zhí)行時(shí)間仍為0．31 s），而是趨近于該額外的時(shí)間開(kāi)銷(xiāo)。因此，在實(shí)際應(yīng)用中，可首先根據(jù)式（9）計(jì)算得到理論上的GPU負(fù)載分配率，然后再減去一個(gè)微小修正量，得到更為精確的分配率。

圖6 GPU對(duì)應(yīng)負(fù)載率變化時(shí)GPU和CPU的執(zhí)行時(shí)間變化情況Fig．6 The variation of GPU and CPU run time with different workload ratio for GPU

對(duì)各版本MTF補(bǔ)償方法執(zhí)行時(shí)間和加速比進(jìn)行統(tǒng)計(jì)，如表3所示。在使用負(fù)載分配策略后，方法執(zhí)行效率進(jìn)一步提高，最終加速比達(dá)到47．82倍，相應(yīng)的執(zhí)行時(shí)間壓縮至1．62 s，可以滿(mǎn)足對(duì)高分辨率衛(wèi)星影像進(jìn)行近實(shí)時(shí)MTF補(bǔ)償?shù)男枨蟆?/p>

表3 各版本MTF補(bǔ)償方法執(zhí)行時(shí)間和加速比Tab．3 The run time and speedup ratio of different MTF compensation approaches

3．2 正確性

為驗(yàn)證本文方法的正確性，將經(jīng)CPU處理得到的影像和由GPU處理得到的影像進(jìn)行逐像素比對(duì)，發(fā)現(xiàn)兩者完全相同，說(shuō)明方法的正確性得到了保證。實(shí)際上，基于Fermi架構(gòu)的GPU的浮點(diǎn)運(yùn)算精度完全滿(mǎn)足IEEE 754－2008浮點(diǎn)運(yùn)算標(biāo)準(zhǔn)［20］，因此對(duì)于需要進(jìn)行浮點(diǎn)運(yùn)算的程序，選擇基于Fermi架構(gòu)的GPU可完全保證計(jì)算結(jié)果的準(zhǔn)確性和精度。

4 結(jié) 論

本文系統(tǒng)地探討了負(fù)載分配的CPU／GPU高分辨率衛(wèi)星影像MTF補(bǔ)償方法。首先闡述了方法的基本實(shí)現(xiàn)，然后通過(guò)3種GPU性能優(yōu)化策略（執(zhí)行配置優(yōu)化、存儲(chǔ)訪問(wèn)優(yōu)化和指令優(yōu)化）對(duì)方法性能進(jìn)行優(yōu)化。使用本文方法在NVIDIA Tesla C2050 GPU上對(duì)高分一號(hào)衛(wèi)星全色影像進(jìn)行MTF補(bǔ)償，與傳統(tǒng)CPU串行方法相比，加速比達(dá)到42．80倍，優(yōu)于文獻(xiàn)［6］中方法。此外，為充分利用CPU的計(jì)算性能，提出了負(fù)載分配策略，該策略根據(jù)系統(tǒng)中計(jì)算核的數(shù)量和方法的實(shí)際加速比為計(jì)算核和GPU分配對(duì)應(yīng)的負(fù)載量，保證方法的總執(zhí)行效率達(dá)到最優(yōu)。使用該策略后，方法執(zhí)行效率進(jìn)一步提高，最終加速比達(dá)到47．82倍，相應(yīng)處理時(shí)間壓縮至1．62 s，可以滿(mǎn)足對(duì)高分辨率衛(wèi)星影像進(jìn)行近實(shí)時(shí)MTF補(bǔ)償?shù)奶幚硇枨?。此外，由于Tesla C2050 GPU的浮點(diǎn)運(yùn)算精度完全滿(mǎn)足IEEE 754－2008浮點(diǎn)運(yùn)算標(biāo)準(zhǔn)，使用本方法處理得到的影像與使用傳統(tǒng)串行CPU方法處理得到的影像完全相同，說(shuō)明了方法的正確性。

［1］ CHEN Qiang，DAI Qiyan，XIA Deshen．Restoration of Remote Sensing Images Based on MTF Theory［J］．Journal of Image and Graphics，2006，11（9）：1299－1305．（陳強(qiáng)，戴奇燕，夏德深．基于MTF理論的遙感影像復(fù)原［J］．中國(guó)圖象圖形學(xué)報(bào)，2006，11（9）：1299－1305．）

［2］ YANG H，DU Q，CHEN G S．Unsupervised Hyperspectral Band Selection Using Graphics Processing Units［J］．IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2011，4（3）：660－668．

［3］ JAVIER R S，MARIA C C，JULIO M H．Geocorrection for Airborne Pushbroom Imagers［J］．IEEE Transactions of Geoscience and Remote Sensing，2012，50（11）：4409－4419．

［4］ MIELIK AINEN J，H U A NG B，H U ANG H L．GPU－accelerated Multi－profile Radiative Transfer Model for the Infrared Atmospheric Sounding Interferometer［J］．IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2011，4（3）：691－699．

［5］ XIAO Han，ZH ANG Zuxun．Parallel Image Matching Algorithm Based on GPGPU［J］．Acta Geodaetica et Cartographica Sinica，2010，39（1）：46－51．（肖漢，張祖勛．基于GPGPU的并行影像匹配算法［J］．測(cè)繪學(xué)報(bào)，2010，39（1）：46－51．）

［6］ FANG L Y，WANG M，LI D R，et al．CPU／GPU Near Realtime Preprocessing for ZY－3 Satellite Images：Relative Radiometric Correction，MTF Compensation and Geocorrection［J］．ISPRS Journal of Photogrammetry and Remote Sensing，2014，87：229－240．

［7］ GE Ping，WANG Mi，PAN Jun，et al．A Study of Adaptive Image Restoration of High Resolution TDI－CCD Image Data［J］．Remote Sensing for Land and Resources，2010（4）：23－28．（葛蘋(píng)，王密，潘俊，等．高分辨率TDI－CCD成像數(shù)據(jù)的自適應(yīng)MTF圖像復(fù)原處理研究［J］．國(guó)土資源遙感，2010（4）：23－28．）

［8］ ZHANG Ying，HE Binbin，LI Xiaowen．Remote Sensing Image Fusion of Beijing－1 DMC＋4 Microsatellite Based on MTF Filter［J］．Acta Geodaetica et Cartographica Sinica，2009，38（3）：223－228．（張瑛，何彬彬，李小文．基于MTF濾波的北京一號(hào)小衛(wèi)星遙感影像融合［J］．測(cè)繪學(xué)報(bào)，2009，38（3）：223－228．）

［9］ LEE D H，YANG J Y，SEO D C，et al．Image Restoration of the Asymmetric Point Spread Function of a Highresolution Remote Sensing Satellite with Time－delayed Integration［J］．Advances in Space Research，2011，47：690－701．

［10］ CAI Xinmin．MTF Calculation and Analysis Based on Satellite Remote Sensing Images［D］．Nanjing：Nanjing University of Science and Technology，2007．（蔡新明．基于衛(wèi)星遙感圖像的 MTF計(jì)算和分析［D］．南京：南京理工大學(xué)，2007．）

［11］ GU Xingfa，LI Xiaoying，MIN Xiangjun，et al．MTF On－orbit Measurement and Image Compensation for CBERS－02 CCD Imager［J］．Science in China Information Sciences，2005，35（1）：26－40．（顧行發(fā)，李小英，閔祥軍，等．CBERS－02衛(wèi)星CCD相機(jī)MTF在軌測(cè)量及圖像MTF補(bǔ)償［J］．中國(guó)科學(xué)E輯：信息科學(xué)，2005，35（1）：26－40．）

［12］ LI Xiaoying，GU Xingfa，YU Tao，et al．Image－derived MTF Method and MTF Compensation for CBERS－02B WFI Imager［J］．Journal of Remote Sensing：2009，13（3）：371－384．（李小英，顧行發(fā)，余濤，等．CBERS－02B衛(wèi)星WFI成像在軌MTF估算與圖像MTF補(bǔ)償［J］．遙感學(xué)報(bào)，2009，13（3）：371－384．）

［13］ KUCK D J，CHAPMAN B，JOST G，et al．Using Open MP：Portable Shared Memory Parallel Programming［M］．Massachusetts：MIT Press，2007．

［14］ YANG Jingyu．Study on Parallel Processing Technologies of Photogrammetry Data Based on GPU［D］．Zhengzhou：Information Engineering University，2011．（楊靖宇．?dāng)z影測(cè)量數(shù)據(jù)GPU并行處理若干關(guān)鍵技術(shù)研究［D］．鄭州：信息工程大學(xué)，2011．）

［15］ ZHANG Shu，CHU Yanli．GPU High Performance Computing—CUDA［M］．Beijing：China Waterpower Press，2009．（張舒，褚艷利．GPU高性能運(yùn)算之CUDA［M］．北京：中國(guó)水利水電出版社，2009．）

［16］ ZHAO Jin．Study of Remote Sensing Image Parallel Processing Algorithms based on GPU and Optimization Techniques［D］．Changsha：The National University of Defense Technology，2011．（趙進(jìn)．基于GPU的遙感圖像并行處理算法及其優(yōu)化技術(shù)研究［D］．長(zhǎng)沙：國(guó)防科學(xué)技術(shù)大學(xué)，2011．）

［17］ LU Fengshun，SONG Junqiang，YIN Fukang et al．Survey of CPU／GPU Synergetic Parallel Computing［J］．Computer Science，2011，38（3）：5－10．（盧風(fēng)順，宋君強(qiáng)，銀?？?，等．CPU／GPU協(xié)同并行計(jì)算技術(shù)研究綜述［J］．計(jì)算機(jī)科學(xué)，2011，38（3）：5－10．）

［18］ KIRK D B，HWU W M W．Programming Massively Parallel Processors［M］．Amsterdam：Elsevier，2010．

［19］ NVIDIA．NVIDIA’s White Paper of Next Generation CUDA Compute Architecture：Fermi［R］．California：NVIDIA Corporation，2009．

［20］ NVIDIA．CUDA：CUFFT Library［S］．California：NVIDIA Corporation，2011．

（責(zé)任編輯：宋啟凡）

A Workload－distribution Based CPU／GPU MTF Compensation Approach for High Resolution Satellite Images

FANG Liuyang，WANG Mi，LI Deren，PAN Jun
State Key Laboratory of Information Engineering in Surveying，Mapping and Remote Sensing，Wuhan University，Wuhan 430079，China

A novel workload－distribution based CPU／GPU MTF compensation approach for high resolution satellite images is proposed．First，the basic GPU implementation issues are addressed；next，three performance tuning methods，the execution configuration optimization，memory access optimization and instruction optimization，are applied to further improve the performance．The approach is tested with the GF－1 panchromatic image in the CPU／GPU system that consists of an Intel Xeon E5650 CPU and a NVIDIA Tesla C2050 GPU．The experimental result shows that the speedup ratio is up to 42．80 times．Furthermore，the CPU／GPU workload distribution strategy is presented to fully exploit CPU’s computing horsepower．With this strategy，the speedup ratio of MTF compensation finally reaches to 47．82 times（the corresponding processing time is 1．62 s），which could meet the requirement of near real－time MTF compensation for high resolution satellite images．

CPU／GPU；MTF compensation；execution configuration optimization；memory access optimization；instruction optimization；workload distribution strategy

FANG Liuyang（1987—），male，PhD candidate，majors in high performance computing in remote sensing．

WANG Mi

P237

A

1001－1595（2014）06－0598－09

國(guó)家973計(jì)劃（2014CB744201；2012CB719902；2012CB719901）；國(guó)家863計(jì)劃（2011AA120203）；中央高?；究蒲袠I(yè)務(wù)費(fèi)專(zhuān)項(xiàng)資金（2012619020205）；長(zhǎng)江學(xué)者與創(chuàng)新團(tuán)隊(duì)發(fā)展計(jì)劃（IRT1278）；國(guó)家自然科學(xué)基金（41371430）

2013－01－29

方留楊（1987—），男，博士生，研究方向?yàn)檫b感數(shù)據(jù)高性能處理。

王密

FANG Liuyang，WANG Mi，LI Deren，et al．A Workload－distribution Based CPU／GPU MTF Compensation Approach for High Resolution Satellite Images［J］．Acta Geodaetica et Cartographica Sinica，2014，43（6）：598－606．（方留楊，王密，李德仁，等．負(fù)載分配的CPU／GPU高分辨率衛(wèi)星影像調(diào)制傳遞補(bǔ)償方法［J］．測(cè)繪學(xué)報(bào)，2014，43（6）：598－606．）

10．13485／j．cnki．11－2089．2014．0095

修回日期：2013－12－13

E－mail：fangliuyang＠whu．edu．cn．

E－mail：wangmi＠whu．edu．cn