趙曉晨,吳皓楠,李林宜,孟令奎
(1.武漢大學(xué)遙感信息工程學(xué)院,武漢 430079;2.珠江水利委員會(huì)珠江水利綜合技術(shù)中心,廣州 510611)
在汛旱情監(jiān)測領(lǐng)域,遙感技術(shù)以其長期動(dòng)態(tài)監(jiān)測的優(yōu)勢,成為強(qiáng)有力的技術(shù)手段之一。遙感對地觀測技術(shù)飛速發(fā)展的同時(shí),帶來了文件體積過大、處理耗時(shí)過長的問題,因此,提高數(shù)據(jù)處理效率、減少耗時(shí)成為亟待解決的問題之一。
汛情及洪澇監(jiān)測應(yīng)用的關(guān)鍵在于水體范圍及水體變化的信息提取,目前常用歸一化差異水體指數(shù)(normalized difference water index,NDWI)、改進(jìn)的歸一化差異水體指數(shù)(modified NDWI,MNDWI)、增強(qiáng)型水體指數(shù)(enhanced water index,WI)等方法提取水體信息,取得了一定成效。旱情監(jiān)測的主要方法之一是利用植被指數(shù)研究不同時(shí)期作物長勢,經(jīng)典的指數(shù)包括歸一化植被指數(shù)(normalized difference vegetation index,NDVI)、比值植被指數(shù)(ratio vegetation index,RVI)、植被狀態(tài)指數(shù)(vegetation condition index,VCI)等。
統(tǒng)一計(jì)算架構(gòu)(compute unified device architecture,CUDA)的基本思想是驅(qū)動(dòng)程序?qū)⒚芗?jì)算數(shù)據(jù)劃分成細(xì)小的區(qū)塊,映射、加載到圖形處理器(graphics processing unit,GPU)中并行執(zhí)行,盡量開發(fā)線程級(jí)并行并將這些線程在硬件中動(dòng)態(tài)調(diào)度和執(zhí)行[1-2]。CUDA采用主機(jī)程序嵌套內(nèi)核函數(shù)的程序設(shè)計(jì)模式,主機(jī)/設(shè)備代碼分離的編譯模式,數(shù)據(jù)由主機(jī)端傳輸至設(shè)備端,并調(diào)用設(shè)備端的內(nèi)核函數(shù)派生大量獨(dú)立線程并行執(zhí)行。GPU與中央處理器(central processing unit,CPU)協(xié)同處理技術(shù)是將GPU視為高性能并行計(jì)算設(shè)備,將CPU稱為主機(jī)端,把GPU看作CPU的協(xié)處理器,避免了從通用算法到應(yīng)用程序接口(application programming interface,API)的復(fù)雜映射,從而直接使用GPU的計(jì)算能力,實(shí)現(xiàn)了異構(gòu)環(huán)境下對GPU資源的靈活調(diào)用。
本文在分析了遙感影像處理原理的基礎(chǔ)上,將處理過程劃分為若干模塊,使用CUDA技術(shù)對影像處理方法做并行改進(jìn),提出了基于CUDA架構(gòu)的CPU-GPU協(xié)同的遙感影像預(yù)處理及相關(guān)指數(shù)計(jì)算并行算法,并選取了典型影像數(shù)據(jù)對該算法進(jìn)行驗(yàn)證,證明其有效性。
汛旱情監(jiān)測應(yīng)用中的遙感影像處理主要包括輻射校正、幾何糾正、遙感指數(shù)計(jì)算等過程,這里對傳統(tǒng)算法原理進(jìn)行分解分析,并確定可優(yōu)化模塊。
輻射校正的目的是將影像像元量化的、無量綱的像元灰度值定標(biāo)為輻射亮度值或表觀反射率,通常定標(biāo)為表觀反射率的計(jì)算表達(dá)式為:
L=GainDN+Bias,
(1)
式中:L為大氣表觀反射率;DN為像元灰度值;Gain為輻射增益參數(shù);Bias為輻射偏置參數(shù),Gain和Bias統(tǒng)稱為輻射定標(biāo)參數(shù)。
實(shí)現(xiàn)輻射定標(biāo)首先要讀取影像頭文件,獲取輻射增益、輻射偏置、太陽高度角等元數(shù)據(jù)信息,其次根據(jù)式(1)對各波段影像柵格進(jìn)行波段運(yùn)算,考慮對波段運(yùn)算處理過程進(jìn)行并行優(yōu)化。
影像幾何糾正的目的是消除圖像的幾何畸變,該過程是遙感數(shù)據(jù)預(yù)處理中耗時(shí)最長、計(jì)算量最大的過程之一,也成為遙感影像快速完成處理的瓶頸。這里采用有理多項(xiàng)式系數(shù)(rational polynomial coefficient,RPC)模型對影像糾正(圖1),相關(guān)研究表明[3],RPC模型相對于嚴(yán)格成像模型的擬合誤差不超過0.04個(gè)像元,均方差小于0.01個(gè)像元,因此可以看作對嚴(yán)格成像模型的擬合。
圖1 RPC模型幾何糾正算法Fig.1 Geometric correction algorithm of RPC model
由圖1知,RPC糾正算法的實(shí)際開發(fā)內(nèi)容主要包括以下4個(gè)模塊:①影像讀取及批處理任務(wù)管理;②RPC正解函數(shù)構(gòu)建;③RPC反解函數(shù)構(gòu)建;④糾正后影像柵格重采樣。
前3個(gè)模塊運(yùn)行效率主要依賴于計(jì)算機(jī)存儲(chǔ)設(shè)備性能,重采樣是RPC糾正中耗時(shí)最長的步驟。使用傳統(tǒng)的串行程序設(shè)計(jì)方式處理時(shí),對所有像素依次循環(huán)處理,造成了極大的流程阻塞,成為算法的最大瓶頸。同時(shí)由于每個(gè)像素的重采樣獨(dú)立進(jìn)行,算法具有可并行性,因此考慮作為主要并行化模塊進(jìn)行設(shè)計(jì)。
1.3.1 汛情遙感監(jiān)測算法
水體范圍變化是汛情及洪澇災(zāi)害的主要特征,水體信息的有效提取也一直是水利遙感技術(shù)研究的重點(diǎn)。常用的水體光譜指數(shù)有NDWI[4]和MNDWI[5]等,具體表達(dá)式分別為:
NDWI=(RGreen-RNIR)/(RGreen+RNIR),
(2)
MNDWI=(RGreen-RMIR)/(RGreen+RMIR),
(3)
式中,RNIR,RGreen和RMIR分別為近紅外、綠光和中紅外波段反射率。
1.3.2 旱情遙感監(jiān)測算法
NDVI,VCI和EVI等遙感指數(shù)可用于監(jiān)測植被長勢,進(jìn)而表征干旱程度。具體表達(dá)式分別為:
NDVI=(RNIR-RRed)/(RNIR+RRed),
(4)
VCI=(NDVI-NDVImin)/(NDVImax-NDVImin),
(5)
(6)
式中:RRed和RBlue分別為紅光和藍(lán)光波段反射率;L為土壤調(diào)節(jié)參數(shù),L=1;C1為紅光大氣改正參數(shù),C1=6;C2為藍(lán)光大氣改正參數(shù),C2=7.5;G為增益系數(shù)。
遙感指數(shù)計(jì)算在算法上無需直方圖統(tǒng)計(jì)、柵格索引變更等復(fù)雜過程,僅包含不同波段影像柵格數(shù)據(jù)中同名像元的波段計(jì)算,這一過程內(nèi)在并行性強(qiáng),因此將波段運(yùn)算模塊作為并行化改造的首要任務(wù)。
遙感影像處理算法并行模式如圖2所示,遙感影像預(yù)處理算法通常以影像數(shù)據(jù)的管理與解析作為起始功能模塊,將影像元數(shù)據(jù)和光譜數(shù)據(jù)柵格分別應(yīng)用于不同的算法模塊。其中元數(shù)據(jù)主要用于相關(guān)模型的構(gòu)建及參數(shù)求解,而光譜數(shù)據(jù)作為影像信息的主要承載,是包括波段運(yùn)算、柵格重采樣、統(tǒng)計(jì)學(xué)分析、矢量化等多種算法模塊的數(shù)據(jù)來源,其中波段運(yùn)算和柵格重采樣具有較強(qiáng)的區(qū)域并行性,適合進(jìn)行并行化改造。同時(shí)不同功能模塊間存在復(fù)雜的協(xié)同關(guān)系,從而使算法理論能夠準(zhǔn)確地應(yīng)用于光譜數(shù)據(jù)本身最終完成產(chǎn)品的輸出。
圖2 遙感影像處理算法并行模式Fig.2 Parallel mode of remote sensing image processing algorithm
根據(jù)1.1節(jié)的原理分析,輻射校正算法的效率瓶頸集中在波段運(yùn)算環(huán)節(jié)。相比于遙感指數(shù)計(jì)算,輻射定標(biāo)涉及的波段更多,公式更復(fù)雜,分支結(jié)構(gòu)更多,計(jì)算密集度更高。從影像頭文件中獲取輻射偏置、輻射增益、太陽高度角等定標(biāo)參數(shù),也是算法的重要環(huán)節(jié)。結(jié)合GPU存儲(chǔ)方式的分析以及定標(biāo)參數(shù)在算法中定位,實(shí)際開發(fā)過程中使用__constant__將參數(shù)映射至常量內(nèi)存,提高算法效率。
首先采用GDAL(geospatial data abstraction library)庫用于數(shù)據(jù)解析,并通過注冊影像驅(qū)動(dòng)的方式,實(shí)現(xiàn)Arc/Info ASCII Grid(asc),GeoTiff (tiff),Erdas Imagine Images(img)和ASCII DEM(dem)等不同格式的影像文件讀取。
影像數(shù)據(jù)解析為柵格數(shù)組后,以數(shù)據(jù)集和波段2級(jí)數(shù)據(jù)類型管理,針對異構(gòu)環(huán)境下設(shè)備端顯存空間有限的問題,采用條帶劃分、分批次處理模式。這種單一的波段運(yùn)算不涉及空間運(yùn)算,僅針對同一像元位置不同波段灰度值,因此對柵格數(shù)組邏輯分塊可以避免多余的空間信息賦值工作,并且根據(jù)GPU的需求彈性調(diào)整分塊大小。
使用核函數(shù)編寫波段計(jì)算是實(shí)現(xiàn)并行處理的核心步驟,實(shí)際編寫中有以下關(guān)鍵點(diǎn):①使用__global__函數(shù)限定符,主機(jī)端負(fù)責(zé)調(diào)用,設(shè)備端負(fù)責(zé)執(zhí)行;②在核函數(shù)中完成柵格數(shù)組索引與線程邏輯管理索引的正確映射;③算法公式中涉及的變量選取合適的變量類型確保分配足夠的顯存空間;③在正確表述算法公式外,使用分支結(jié)構(gòu)處理函數(shù)異常值。
其中線程管理索引決定了柵格數(shù)組元素在線程邏輯結(jié)構(gòu)中的分配方式,其公式為:
I=(byGx+bx)(BxBy)+tyBx+tx,
(7)
式中:I為柵格數(shù)組索引;bx和by分別為線程塊在網(wǎng)格中x和y方向上的索引,對應(yīng)內(nèi)建變量blockIdx.x和blockIdx.y;Gx為一個(gè)網(wǎng)格中x方向上線程塊總數(shù),對應(yīng)內(nèi)建變量gridDim.x;Bx和By分別為一個(gè)線程塊中x和y方向上線程總數(shù),對應(yīng)內(nèi)建變量blockDim.x和blockDim.y;tx和ty分別為線程在線程區(qū)塊中x和y方向上的索引,對應(yīng)內(nèi)建變量threadIdx.x和threadIdx.y。本文結(jié)合二維影像數(shù)據(jù)自身特點(diǎn),同樣使用二維格網(wǎng)、線程塊進(jìn)行管理。
另外,由于核函數(shù)中采用了少量分支結(jié)構(gòu)用于波段運(yùn)算中異常值的處理,可能造成線程塊中執(zhí)行速率不一致。為此可使用CUDA中內(nèi)建函數(shù)cudaThreadSynchronize()保證同一線程塊中各線程狀態(tài)對齊。
綜上,針對波段運(yùn)算的并行算法設(shè)計(jì)中,主機(jī)端負(fù)責(zé)批處理任務(wù)編排、影像文件讀取和輸出、影像柵格數(shù)組劃分等耗時(shí)較短的串行任務(wù),設(shè)備端負(fù)責(zé)對影像柵格數(shù)據(jù)并行執(zhí)行核函數(shù),完成指數(shù)算法及異常值處理。這一設(shè)計(jì)也是RPC幾何糾正并行算法的基礎(chǔ)。
通過1.2節(jié)中對基于RPC模型的幾何糾正算法的原理分析,結(jié)合CUDA程序設(shè)計(jì)模型設(shè)計(jì)了CPU-GPU協(xié)同的RPC糾正算法,如圖3所示。
圖3 并行化RPC糾正算法Fig.3 Parallel RPC algorithm
該算法實(shí)施步驟如下:①影像解析及參數(shù)讀取,并存儲(chǔ)于主機(jī)端內(nèi)存;②根據(jù)RPC模型構(gòu)建正反解函數(shù),正解函數(shù)在主機(jī)端調(diào)用并執(zhí)行,且僅用于糾正后影像四至點(diǎn)的結(jié)算,反解函數(shù)將用于設(shè)備端影像并行重采樣;③矯正后影像柵格初始化,填充重采樣灰度數(shù)據(jù),同時(shí)對該范圍進(jìn)行邏輯分塊,以便在有限的設(shè)備端存儲(chǔ)空間下完成整張影像的重采樣,分塊后影像條帶范圍被傳入GPU,在全局內(nèi)存中初始化地面點(diǎn)坐標(biāo)數(shù)組,根據(jù)校正后影像條帶范圍反解得到相對應(yīng)的原始影像范圍,將該范圍內(nèi)原始影像柵格傳入GPU;④RPC反解函數(shù)中各項(xiàng)系數(shù)隨第一個(gè)影像塊傳入GPU,并使用__constant__限定符存儲(chǔ)為常量內(nèi)存,然后將反解函數(shù)派生為多線程任務(wù),將分塊影像條帶范圍內(nèi)的地面點(diǎn)坐標(biāo)并行解算為對應(yīng)的影像坐標(biāo),并根據(jù)此坐標(biāo)對原始影像條帶進(jìn)行重采樣;⑤將糾正后影像條帶回傳至主機(jī)端,合并入糾正后影像柵格。最后使用GDAL庫中的Geo-TIFF影像驅(qū)動(dòng)封裝為影像產(chǎn)品。其中用于并行坐標(biāo)反解及重采樣主要包括并行糾正核函數(shù)、RPC模型填充函數(shù)以及RPC反解函數(shù)等幾類。
遙感指數(shù)并行算法與輻射校正并行算法相似,對讀取后的影像波段直接進(jìn)行數(shù)學(xué)計(jì)算。對于多波段影像,相較于輻射定標(biāo)過程,指數(shù)計(jì)算涉及的波段數(shù)更少,因此也更為簡單。
在2.1節(jié)波段運(yùn)算并行算法設(shè)計(jì)中提到了對柵格數(shù)組進(jìn)行邏輯劃分以避免顯存溢出,以及柵格數(shù)組在多級(jí)線程管理結(jié)構(gòu)中的映射方法,這2項(xiàng)算法設(shè)計(jì)實(shí)現(xiàn)了異構(gòu)環(huán)境下高空間分辨率影像分批并行處理,是提升并行算法效率的關(guān)鍵。因此這里結(jié)合GPU硬件參數(shù)和相關(guān)官方資料進(jìn)行了實(shí)驗(yàn),并對相關(guān)算法參數(shù)進(jìn)行優(yōu)化。
線程束(warp)作為GPU執(zhí)行任務(wù)時(shí)的調(diào)度單位,由調(diào)度器分配依次進(jìn)入流處理器簇(SM)執(zhí)行。SM每次只執(zhí)行線程塊(block)中的一個(gè)warp,稱作active warp,此時(shí)其余warp處于待命狀態(tài)。當(dāng)active warp在執(zhí)行中需要等待時(shí)(如訪問GPU全局內(nèi)存),調(diào)度器可以直接切換至其他warp執(zhí)行。通過這種方式GPU隱藏了多線程的延遲等待,實(shí)現(xiàn)了大規(guī)模的并行。因此在程序設(shè)計(jì)時(shí),可以在同一block設(shè)置盡量多的warp,避免所有warp均處于等待訪存狀態(tài)。但由于SM內(nèi)寄存器數(shù)量有限,同一個(gè)block內(nèi)warp數(shù)量過多可能造成單一線程可使用的寄存器數(shù)量過少,影響數(shù)據(jù)訪存效率。同時(shí)由于warp中線程(thread)數(shù)為固定值,因此同一block內(nèi)thread總數(shù)應(yīng)設(shè)計(jì)為這一固定值的倍數(shù),避免warp切換時(shí)流處理器(SP)的閑置。通過查閱相關(guān)文檔獲得了實(shí)驗(yàn)平臺(tái)所使用GPU相關(guān)硬件參數(shù)如表1所示。
表1 GPU硬件參數(shù)Tab.1 GPU parameters
GPU硬件及CUDA詳細(xì)參數(shù)如表2。實(shí)驗(yàn)使用便攜計(jì)算機(jī)平臺(tái),CPU版本為Intel(R)Core(TM)i5-2520M CPU @ 2.50 GHz,GPU為NVIDIA GeForce GT 640 M LE,DDR3內(nèi)存為16 GB。影像數(shù)據(jù)Landsat8 Level-1產(chǎn)品,OLI和TIRS傳感器,共11個(gè)波段,覆蓋周期為16 d,掃描寬度為170 km×180 km,單波段影像數(shù)據(jù)量為100 MB。高分一號(hào)(GF-1)L1A產(chǎn)品,共9個(gè)波段,PMS傳感器幅寬為60 km,空間分辨率為2~8 m,WFV傳感器幅寬為800 km,空間分辨率為16m,單景影像數(shù)據(jù)量為1.2 GB。
表2 GPU硬件及CUDA詳細(xì)參數(shù)Tab.2 GPU and CUDA parameters
根據(jù)3.1節(jié)中的硬件參數(shù),在控制實(shí)驗(yàn)平臺(tái)其他硬件參數(shù)不變的前提下,使用1景GF-1 WFV影像(數(shù)據(jù)量1.63 GB)測試了不同線程塊尺寸(以32為基準(zhǔn))下NDVI并行處理算法中波段運(yùn)算環(huán)節(jié)處理速度,結(jié)果如圖4所示。
圖4 線程塊尺寸實(shí)驗(yàn)Fig.4 Different thread block sizes
可見在線程塊內(nèi)最大線程數(shù)范圍內(nèi),線程塊尺寸對波段運(yùn)算耗時(shí)的影響并不顯著。原因可能在于本算法涉及的數(shù)據(jù)類型較為單一,各線程對全局內(nèi)存訪存延遲較短。
對于可能影響算法運(yùn)行效率的另一環(huán)節(jié)——柵格數(shù)組邏輯分塊,同樣控制了其他變量不變,以影像寬度為基準(zhǔn),分別測試了不同行數(shù)的條帶劃分尺度下波段運(yùn)算耗時(shí),實(shí)驗(yàn)結(jié)果如圖5所示。
圖5 柵格劃分尺度實(shí)驗(yàn)Fig.5 Grid scale of performance
由圖5可以看出,柵格數(shù)組的劃分尺度的確影響了并行波段運(yùn)算的效率,更大的條帶劃分可以減少CPU-GPU間數(shù)據(jù)交換的頻數(shù),縮短并行算法的時(shí)間。但條帶劃分仍需考慮GPU顯存容量的上限,結(jié)合實(shí)驗(yàn)和理論分析可以得出最佳條帶劃分公式為:
(8)
式中:T為最佳條帶寬度;H和W分別為二維影像柵格數(shù)組的高度和寬度;n為算法使用的柵格數(shù)組數(shù)量;M為GPU全局內(nèi)存;m為影像波段數(shù)量。
3.3.1 汛旱情監(jiān)測指數(shù)
分別使用基于CUDA的并行算法、傳統(tǒng)串行算法處理一景GF-1 WFV影像數(shù)據(jù),進(jìn)行NDVI和NDWI指數(shù)產(chǎn)品的生產(chǎn),對各模塊計(jì)算耗時(shí)統(tǒng)計(jì)后分別見圖6。以鄱陽湖為例計(jì)算光譜指數(shù)產(chǎn)品如圖7所示。
(a)NDVI指數(shù)算法耗時(shí)對比 (b)NDWI指數(shù)算法耗時(shí)對比
(a)NDVI (b)NDWI
對于簡單波段運(yùn)算類光譜指數(shù)算法,采用基于CUDA的并行算法能夠有效提高波段運(yùn)算這一瓶頸環(huán)節(jié)的效率,影像讀取和產(chǎn)品輸出環(huán)節(jié)的時(shí)間消耗與串行算法基本一致。此外由于波段運(yùn)算核函數(shù)中具體內(nèi)容與串行算法循環(huán)體中函數(shù)內(nèi)容完全一致,因此保證了產(chǎn)品結(jié)果的正確性。
3.3.2 輻射校正
分別使用基于CUDA的并行算法、傳統(tǒng)串行算法對1景Landsat8影像數(shù)據(jù)進(jìn)行輻射定標(biāo)。由于Landsat8各波段影像使用單獨(dú)的TIFF文件存儲(chǔ),因此使用各波段處理的平均耗時(shí)作為對比。結(jié)果如表3所示。表3中可見,由于各波段影像差異和計(jì)算機(jī)存在性能因素,算法耗時(shí)存在一定范圍的浮動(dòng)。從整體來看,波段運(yùn)算環(huán)節(jié)處理效率有了大幅提升,相比于串行處理方法,并行算法可節(jié)約58.9%的時(shí)間,但直方圖統(tǒng)計(jì)環(huán)節(jié)仍然有一定的優(yōu)化空間,算法總耗時(shí)得到了一定程度的降低。
表3 傳統(tǒng)串行和CUDA并行輻射處理耗時(shí)對比Tab.3 Processing time of CPU and CUDA for radiometric correction (s)
(續(xù)表)
3.3.3 RPC幾何糾正
分別使用基于CUDA的并行RPC糾正算法、傳統(tǒng)串行RPC糾正算法,對1景GF-1 WFV影像數(shù)據(jù)進(jìn)行幾何糾正。在影像重采樣環(huán)節(jié)使用了最鄰近像元法和雙線性內(nèi)插法,比較不同重采樣算法對并行加速比的影響,如圖8和表4所示。根據(jù)圖8和表4可以發(fā)現(xiàn),基于CUDA的并行重采樣算法相比傳統(tǒng)串行算法可以實(shí)現(xiàn)7~10倍的加速比,且不同重采樣方式的加速效果較為接近。使用較為復(fù)雜重采樣方式總處理時(shí)間較長、并行處理占比較大,因此總加速比更高。
(a)最鄰近像元重采樣 (b)雙線性內(nèi)插重采樣
表4 RPC并行校正算法重采樣環(huán)節(jié)加速比Tab.4 Speedup radio of parallel RPC algorithm
針對汛旱情監(jiān)測數(shù)據(jù)處理包括輻射校正、幾何糾正及遙感指數(shù)計(jì)算等過程耗時(shí)過長的問題,面向常見遙感數(shù)據(jù)如Landsat和GF-1等,分析各類處理過程的共性與性能瓶頸,基于CUDA架構(gòu)提出了基于功能模塊劃分的并行處理算法。
通過實(shí)驗(yàn)發(fā)現(xiàn),調(diào)整核心參數(shù)進(jìn)行算法訪存優(yōu)化,可以有效提升并行運(yùn)算的效率,實(shí)驗(yàn)結(jié)果表明該算法在保證成果精度的同時(shí),相比傳統(tǒng)串行算法效率有了顯著提升。輻射校正與指數(shù)計(jì)算均涉及的波段計(jì)算模塊,經(jīng)過并行優(yōu)化后,較串行模式可減少70%以上的時(shí)間;幾何糾正不同重采樣方式的加速效果較為接近;通過調(diào)整核心參數(shù)進(jìn)行算法訪存優(yōu)化,可以有效提升并行波段運(yùn)算的效率。本算法針對便攜化處理平臺(tái)而設(shè)計(jì),可為汛旱情監(jiān)測提供切實(shí)有效的高性能計(jì)算解決方案。