陸冬華，趙英俊，張東輝，秦 凱

（1. 核工業(yè)北京地質(zhì)研究院 遙感信息與圖像分析技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室，北京 100029）

遙感對(duì)地觀測(cè)作為一項(xiàng)重要的空間探測(cè)技術(shù)，自20 世紀(jì)60 年代以來逐漸發(fā)展，已在氣象、農(nóng)業(yè)、林業(yè)、資源、環(huán)境、生態(tài)、水利、海洋、大氣、水文、災(zāi)害、全球變化等民用領(lǐng)域，以及情報(bào)偵查、戰(zhàn)場(chǎng)監(jiān)視、導(dǎo)彈制導(dǎo)、隱身等軍用領(lǐng)域發(fā)揮了積極作用[1]。遙感傳感器技術(shù)的發(fā)展，使得獲取的遙感數(shù)據(jù)空間分辨率和光譜分辨率不斷提高，數(shù)據(jù)量不斷增大，因此需要更強(qiáng)的計(jì)算資源對(duì)其進(jìn)行處理[2]。美國(guó)最先進(jìn)的WorldView3衛(wèi)星遙感數(shù)據(jù)，其全色空間分辨率為0.3 m，多光譜16個(gè)波段的空間分辨率為1.24 m[3]，每100 km2融合后的數(shù)據(jù)約為17 GB；加拿大Itres公司的航空高光譜傳感器CASI和SASI，分別擁有32和100個(gè)波段[4]，相對(duì)飛行高度為1 000 m時(shí)，空間分辨率分別為0.49 m和1.22 m，長(zhǎng)度為10 km的單條測(cè)線獲取的數(shù)據(jù)分別約為8 GB和15 GB。目前大多數(shù)的商業(yè)軟件均未采用多線程數(shù)據(jù)處理技術(shù)，因此對(duì)數(shù)據(jù)量巨大的遙感圖像進(jìn)行分析時(shí)需要耗費(fèi)大量的計(jì)算時(shí)間。

國(guó)內(nèi)外學(xué)者對(duì)多線程和并行處理遙感圖像方面開展了相關(guān)研究工作。QIN C Z[5]等利用多線程的方式對(duì)存儲(chǔ)在文件系統(tǒng)中的遙感數(shù)據(jù)進(jìn)行并行讀寫，降低了數(shù)據(jù)讀寫在遙感影像處理中占用的時(shí)間；劉世永[6]等利用MPI技術(shù)在超微高性能計(jì)算機(jī)上實(shí)現(xiàn)了多線程原始柵格影像的標(biāo)準(zhǔn)瓦片輸出；吳煒[7]等針對(duì)集群計(jì)算機(jī)實(shí)現(xiàn)了高分辨率遙感影像的多尺度分割；周海芳[8]等實(shí)現(xiàn)了基于GPU的圖像配準(zhǔn)并行計(jì)算；LYU Z H[9]等針對(duì)高性能集群實(shí)現(xiàn)了遙感圖像聚類分析。

遙感圖像并行處理面臨兩個(gè)難題[5]：遙感數(shù)據(jù)量巨大，從計(jì)算機(jī)硬盤調(diào)入內(nèi)存需耗費(fèi)大量時(shí)間；存儲(chǔ)遙感數(shù)據(jù)的文件格式種類繁多。目前國(guó)內(nèi)外的研究主要側(cè)重于多臺(tái)計(jì)算機(jī)或集群的分布式計(jì)算以及GPU并行計(jì)算，該類方法雖可大幅提高數(shù)據(jù)的處理效率，但由于大多數(shù)科研人員仍以單臺(tái)計(jì)算機(jī)運(yùn)算為主要工作方式，且編程復(fù)雜度高，分布式需綜合考慮不同計(jì)算機(jī)之間的數(shù)據(jù)傳輸與協(xié)同問題，GPU的程序開發(fā)需綜合考慮GPU內(nèi)不同類型顯存的性能差異以及顯存與內(nèi)存之間傳輸數(shù)據(jù)的低效率問題[10]。隨著計(jì)算機(jī)硬件的發(fā)展，Intel 和AMD等多核處理器的市場(chǎng)占有率逐步上升，而以往的分類算法一般為串行執(zhí)行，無法發(fā)揮多核處理器的計(jì)算優(yōu)勢(shì)[11]。因此，本文重點(diǎn)針對(duì)單機(jī)多核環(huán)境下的遙感圖像并行處理關(guān)鍵技術(shù)進(jìn)行研究，一方面利用空間數(shù)據(jù)抽象庫(kù)（GDAL）開源庫(kù)解決不同格式數(shù)據(jù)的輸入輸出問題，另一方面通過OpenMP解決單臺(tái)計(jì)算機(jī)CPU資源的充分利用問題，降低了開發(fā)難度，有效提高了圖像處理速度。

1 不同類型遙感數(shù)據(jù)的輸入和輸出

由于多源遙感數(shù)據(jù)來自不同的衛(wèi)星中心和傳感器，遙感數(shù)據(jù)格式缺乏統(tǒng)一標(biāo)準(zhǔn)，不同的傳感器數(shù)據(jù)具有不同的文件格式，甚至同一傳感器，因其來源不同，文件格式也略有差異[12]。為了解決這一問題，1998年加拿大的Warmerdam F開始了GDAL項(xiàng)目的編寫工作，并得到了許多個(gè)人和團(tuán)體的支持[13]。截至目前GDAL已能支持超過150種空間數(shù)據(jù)格式，并被ESRI的ArcGIS10、Google Earth，跨平臺(tái)的GRASS GIS 系統(tǒng)等多家遙感和GIS軟件使用[14]。

GDAL 是一個(gè)獨(dú)立的專業(yè)開源庫(kù)，是在X/MIT 許可協(xié)議下的開源柵格空間數(shù)據(jù)轉(zhuǎn)換庫(kù)[13]。GDAL 的特點(diǎn)為[15]：①可擴(kuò)展性、可移植性較好，若想支持某種新的數(shù)據(jù)格式，只需添加該格式驅(qū)動(dòng)程序即可；②API 功能全面，可嵌入到其他程序中；③可跨平臺(tái)運(yùn)行；④支持常用的空間數(shù)據(jù)格式。

GDAL讀寫柵格數(shù)據(jù)最主要的核心類如圖1所示。GDAL讀寫文件操作流程如圖2所示。

GDALDataset類繼承自GDALMajorObject，負(fù)責(zé)管理存儲(chǔ)在硬盤上不同文件格式的遙感數(shù)據(jù)，利用GDALDataset中的函數(shù)可獲取遙感數(shù)據(jù)的元數(shù)據(jù)信息、投影坐標(biāo)信息、波段數(shù)量等；并可利用該類中的RasterIO函數(shù)完成多波段數(shù)據(jù)讀寫操作。

GDALDriver類是數(shù)據(jù)格式的抽象，每種格式均會(huì)生成GDALDriver的一個(gè)實(shí)例，所有屬性和方法均來自于其對(duì)應(yīng)的GDALDataset類[16]，通常用于創(chuàng)建不同格式的遙感數(shù)據(jù)文件。將遙感數(shù)據(jù)格式的名稱傳遞給GDALDriver類，其將裝載相應(yīng)格式的信息，并創(chuàng)建相應(yīng)的GDALDataset類，實(shí)現(xiàn)對(duì)文件的讀寫。

GDALRasterBand類負(fù)責(zé)管理GDALDataset類中的單一波段，獲取該波段的信息和數(shù)據(jù)，實(shí)現(xiàn)數(shù)據(jù)的讀寫操作。

圖2 GDAL讀寫文件操作流程圖

2 OpenMP共享存儲(chǔ)編程

OpenMP 是由世界上主要的計(jì)算機(jī)硬件和軟件廠商提出的標(biāo)準(zhǔn)，跨平臺(tái)、可伸縮的模型為并行程序的設(shè)計(jì)者在桌面電腦和巨型機(jī)上開發(fā)并行程序提供了簡(jiǎn)單靈活的接口，規(guī)范了一系列編譯指導(dǎo)語(yǔ)句、子程序庫(kù)和環(huán)境變量，在不降低性能的情況下克服了機(jī)器專用編譯指導(dǎo)帶來的移植性問題，得到了各界的認(rèn)可。自發(fā)布以來，OpenMP已得到工業(yè)界和學(xué)術(shù)界的大力推動(dòng)，成為共享主存多處理機(jī)上并行編程的事實(shí)標(biāo)準(zhǔn)。理論上，OpenMP更好地利用了共享內(nèi)存體系結(jié)構(gòu)，允許運(yùn)行時(shí)調(diào)度，避免了消息傳遞的開銷，并提供了細(xì)粒度和粗粒度運(yùn)行機(jī)制[17-18]。

并行計(jì)算的主要實(shí)現(xiàn)過程為：主進(jìn)程（處理器）將一個(gè)大任務(wù)（待處理數(shù)據(jù)）分派到若干個(gè)子進(jìn)程（處理器）上進(jìn)行處理，再由主進(jìn)程負(fù)責(zé)收集不同子進(jìn)程的數(shù)據(jù)處理結(jié)果并進(jìn)行組合，以達(dá)到多處理器共同完成某一任務(wù)的目的[19]。首先由主線程執(zhí)行程序的串行部分，再由派生出的其他線程執(zhí)行程序的并行部分。相對(duì)于其他并行算法而言，OpenMP并行算法的開發(fā)十分簡(jiǎn)單，即便是已開發(fā)的單線程算法也十分容易改造，無需程序員進(jìn)行復(fù)雜的線程創(chuàng)建、同步、負(fù)載平衡和銷毀等工作。

OpenMP多線程的編寫非常簡(jiǎn)單，包括兩種方式：①在循環(huán)體for(int i=0;i＜N;i++)結(jié)構(gòu)的前一行加入#pragma omp parallel for標(biāo)識(shí)；②通過#pragma omp parallel{…} 標(biāo)記并行塊范圍。系統(tǒng)將根據(jù)當(dāng)前線程數(shù)自動(dòng)對(duì)任務(wù)進(jìn)行分配，利用omp_set_num_threads(int num)可設(shè)置并行區(qū)中活動(dòng)的線程個(gè)數(shù)。使用OpenMP時(shí)，需在包含文件中加入函數(shù)庫(kù)＜o(jì)mp.h＞。常用的函數(shù)和編譯指導(dǎo)語(yǔ)句如表1[20-21]所示。

表1 OpenMP常用函數(shù)和編譯指導(dǎo)語(yǔ)句

3 遙感圖像單機(jī)并行化處理

遙感圖像并行化處理的前提是遙感圖像各部分計(jì)算不存在依賴關(guān)系，即在串行計(jì)算時(shí)，先計(jì)算得到的結(jié)果不對(duì)后計(jì)算得到的結(jié)果產(chǎn)生影響。遙感所涉及的算法種類繁多，大多可實(shí)現(xiàn)并行化處理，如圖像變換、圖像配準(zhǔn)、圖像融合、目標(biāo)提取等；但各方法處理過程完全不同，因此并行化算法的開發(fā)也不盡相同。為了降低算法的并行化難度，突出遙感并行化算法的一般過程，本文以Majority濾波和局部均方差計(jì)算為例進(jìn)行闡述。

3.1 遙感圖像的分塊

由于遙感圖像輸入輸出數(shù)據(jù)占用了大量的磁盤空間，無法一次性將所有數(shù)據(jù)讀入內(nèi)存，因此必須將其進(jìn)行分塊處理，每次僅讀取數(shù)據(jù)的一個(gè)分塊到內(nèi)存中，處理后再輸出到磁盤空間。分塊方法通常包括行優(yōu)先、列優(yōu)先和平均分塊[5]3種。行優(yōu)先處理方法，根據(jù)內(nèi)存大小，將遙感圖像分成行數(shù)相近的若干個(gè)數(shù)據(jù)塊；列優(yōu)先處理方法，根據(jù)內(nèi)存大小，將遙感圖像分成列數(shù)相近的若干個(gè)數(shù)據(jù)塊；平均分塊處理方法，根據(jù)內(nèi)存大小，將遙感圖像分成大小相近的數(shù)據(jù)塊，每個(gè)數(shù)據(jù)塊中數(shù)據(jù)的行數(shù)和列數(shù)相當(dāng)。

處理過程中，由于Majority濾波和局部均方差計(jì)算均需對(duì)M×M窗口內(nèi)的像素值進(jìn)行統(tǒng)計(jì)和計(jì)算，為了保證分塊數(shù)據(jù)邊界計(jì)算的正確，輸入的數(shù)據(jù)應(yīng)在分塊的基礎(chǔ)上，在上下左右各增加(int)(M/2)，即若第i個(gè)分塊的左上角坐標(biāo)為(xoff, yoff)，大小為(xBlockSize,yBlockSize)，則讀取的數(shù)據(jù)范圍應(yīng)為(xoff-(int)(M/2),yoff-(int)(M/2))，大小應(yīng)為（xBlockSize+M－1,yBlockSize+M－1），其中M為基數(shù)。

3.2 遙感并行算法的實(shí)現(xiàn)

首先利用GDAL打開待處理的圖像文件，讀取圖像的行列數(shù)、坐標(biāo)信息、投影信息等基本信息；再計(jì)算分塊大小，確定每個(gè)分塊的左上角坐標(biāo)以及行數(shù)和列數(shù)，利用OpenMP開始多線程計(jì)算，利用GDAL讀取分塊數(shù)據(jù)，并進(jìn)行Majority計(jì)算；然后利用GDAL寫入處理后的數(shù)據(jù)；當(dāng)所有線程都結(jié)束后，關(guān)閉遙感數(shù)據(jù)文件。

以數(shù)據(jù)分塊結(jié)果構(gòu)建for循環(huán)體，并在for循環(huán)體前加入#pragma omp parallel for編譯指導(dǎo)語(yǔ)句，由于在每個(gè)線程均存在IO操作，而單機(jī)多線程的讀寫將顯著降低讀寫效率，因此在IO語(yǔ)句前加入#pragma omp critical原子處理編譯指導(dǎo)語(yǔ)句，避免多線程同時(shí)進(jìn)行IO操作。for循環(huán)體前的編譯指導(dǎo)語(yǔ)句為：

//#pragma omp parallel for

for (int iBlock= 0; iBlock ＜ 分塊總數(shù) ;++ iBlock)

IO操作前的編譯指導(dǎo)語(yǔ)句為：

#pragma omp critical //原子處理

iBand-＞RasterIO(GF_Write, xOff, yOff, xBlockSize,yBlockSize,……);

4 數(shù)據(jù)實(shí)驗(yàn)與分析

4.1 實(shí)驗(yàn)算法

1）Majority算法，主要應(yīng)用于遙感圖像分類的后處理。利用M×N的移動(dòng)窗口，對(duì)窗口內(nèi)的圖像像素值進(jìn)行統(tǒng)計(jì)，將出現(xiàn)頻率最多的像素值設(shè)置為M×N窗口的中心像素值，M和N通常為基數(shù)。

2）局部均方差算法，主要應(yīng)用于圖像增強(qiáng)和邊緣檢測(cè)。利用M×N的移動(dòng)窗口，對(duì)窗口內(nèi)的圖像像素值進(jìn)行標(biāo)準(zhǔn)方差計(jì)算。

4.2 數(shù)據(jù)實(shí)驗(yàn)

實(shí)驗(yàn)平臺(tái)為惠普計(jì)算機(jī)圖形工作站Z800，內(nèi)存為64 GB，2顆Intel Xeon x5690 CPU，主頻為3.46 GHz，每顆CPU為6核，開發(fā)平臺(tái)為Viusal Studio，將項(xiàng)目→屬性→VC++目錄→包含目錄和庫(kù)目錄中分別設(shè)置為GDAL的頭文件目錄和LIB文件目錄。將GDAL的Bin文件夾設(shè)置于環(huán)境變量中，或直接拷貝Bin文件夾中的文件到工程數(shù)據(jù)文件目錄，完成對(duì)GDAL庫(kù)編譯環(huán)境的設(shè)置；然后在C/C++的語(yǔ)言選項(xiàng)中設(shè)置OpenMP支持為“是(/openmp)”，完成對(duì)OpenMP編譯環(huán)境的設(shè)置。

以Majority算法為實(shí)驗(yàn)對(duì)象，輸入圖像大小為3 000×2 500的單波段分類數(shù)據(jù)，輸出為Majority濾波結(jié)果，分別測(cè)試在窗口尺寸為3×3和5×5，并發(fā)線程數(shù)為1～12情況下的計(jì)算性能，如表2所示，可以看出，線程并發(fā)數(shù)由1增加到4的過程中，計(jì)算時(shí)間逐漸縮短，當(dāng)并發(fā)線程增加到5之后，計(jì)算時(shí)間開始變長(zhǎng)，達(dá)到12時(shí)計(jì)算時(shí)間最長(zhǎng)。圖3更加清晰地反映了這一過程。

表2 Majority算法不同線程耗時(shí)表

圖3 Majority算法不同線程耗時(shí)曲線

Majority算法實(shí)驗(yàn)說明盡管計(jì)算機(jī)CPU共有12個(gè)核心，但計(jì)算時(shí)間并未隨并發(fā)數(shù)量的增加而逐漸縮短，而是在4個(gè)線程并發(fā)時(shí)達(dá)到最高效率。為了進(jìn)一步研究并行算法的效率，本文繼續(xù)對(duì)局部均方差算法進(jìn)行實(shí)驗(yàn)。

以局部均方差算法為實(shí)驗(yàn)對(duì)象，輸入圖像大小為3 000×2 500的單波段16 bit整型灰度數(shù)據(jù)，輸出為方差計(jì)算結(jié)果，分別測(cè)試在窗口尺寸為3×3和11×11，并發(fā)線程數(shù)為1～12情況下的計(jì)算性能，具體如表3所示，可以看出，從線程并發(fā)數(shù)由1增加到11的過程中，計(jì)算時(shí)間基本上逐漸縮短，達(dá)到12后略有延長(zhǎng)。圖4更加清晰地反映了這一過程。

表3 局部均方差算法不同線程耗時(shí)表

圖4 局部均方差算法不同線程耗時(shí)曲線

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

實(shí)驗(yàn)結(jié)果表明，局部均方差算法通過OpenMP的加速效果較好；而當(dāng)并發(fā)線程數(shù)為4時(shí)，Majority算法的計(jì)算效率最高。為了進(jìn)一步分析導(dǎo)致實(shí)驗(yàn)結(jié)果的原因，對(duì)參與計(jì)算的各函數(shù)進(jìn)行逐一計(jì)算，分析各函數(shù)對(duì)并行效率的影響。其代碼主要結(jié)構(gòu)為：

omp_set_num_threads(nThread);

#pragma omp parallel for

for(int i= 0;i＜M;i++)

for(int j=0;j＜N;j++)

{

//測(cè)試代碼

…

}

當(dāng)測(cè)試代碼為空或僅包含基本的數(shù)學(xué)運(yùn)算（加、減、乘、除、開方運(yùn)算）時(shí)，并發(fā)線程由1增加至12時(shí)，程序運(yùn)行時(shí)間逐漸縮短。

當(dāng)測(cè)試代碼為int a[5000]或int *a=new int[5000/nThread]時(shí)，并發(fā)線程由1增加至6時(shí)，程序運(yùn)行時(shí)間逐漸縮短，增加至7時(shí)，運(yùn)行時(shí)間有所延長(zhǎng)，但隨后又逐漸縮短。

當(dāng)測(cè)試代碼為int *a= new int[10]時(shí)，程序運(yùn)行時(shí)間與并發(fā)線程數(shù)的變化與Majority算法的規(guī)律相同，即1～4個(gè)線程，程序運(yùn)行時(shí)間逐漸縮短，之后隨并發(fā)數(shù)增加程序運(yùn)行時(shí)間逐漸延長(zhǎng)。

Majority算法中使用了C++STL map類，該類對(duì)并發(fā)線程的加速也有較大影響，當(dāng)測(cè)試代碼部分為map＜int,int＞ statics變量聲明時(shí)，1～5個(gè)并發(fā)線程數(shù)時(shí)，程序運(yùn)行時(shí)間縮短，隨后延長(zhǎng)；將map＜int,int＞ statics變量聲明放在for循環(huán)體之外，測(cè)試代碼僅為statics.find(2)時(shí)，1～4個(gè)線程程序運(yùn)行時(shí)間逐漸縮短，之后隨并發(fā)數(shù)增加程序運(yùn)行時(shí)間逐漸延長(zhǎng)。

由此可知，影響Majority算法加速效果的主要原因?yàn)镃++STL map類的變量聲明和查找函數(shù)；由map類中的哪個(gè)函數(shù)產(chǎn)生的影響還需進(jìn)一步研究，但比較明確的是，程序在運(yùn)行過程中申請(qǐng)內(nèi)存或新建數(shù)組，會(huì)在一定程度上影響OpenMP的加速效果。

5 結(jié) 語(yǔ)

本文將GDAL開源庫(kù)與OpenMP相結(jié)合，開展了遙感數(shù)據(jù)在單機(jī)多核計(jì)算機(jī)環(huán)境下的處理實(shí)驗(yàn)，利用GDAL實(shí)現(xiàn)了多種遙感數(shù)據(jù)的讀寫操作，利用OpenMP將串行程序改寫為并行程序。GDAL讀取遙感數(shù)據(jù)后，應(yīng)根據(jù)內(nèi)存大小對(duì)遙感數(shù)據(jù)進(jìn)行分塊處理，利用OpenMP將for循環(huán)體并行化，可有效提高遙感圖像的處理效率。

實(shí)驗(yàn)結(jié)果表明，并行程序?qū)b感圖像處理的加速效果與處理算法緊密相關(guān)，不同的算法需要調(diào)整并發(fā)線程數(shù)才能達(dá)到最佳的處理效率，單純的增加線程有時(shí)不但不能提高處理效率，反而會(huì)延長(zhǎng)處理時(shí)間。