衛(wèi)星SAR數(shù)據(jù)讀取與保存方法研究與軟件實現(xiàn)

劉元廷，趙朝方，王志雄，馬佑軍

（中國海洋大學(xué) 海洋遙感研究所，山東 青島 266100）

合成孔徑雷達（Synthetic Aperture Radar，SAR）作為主動微波傳感器可以全天時、全天候的工作，具有良好的空間分辨率，因此在環(huán)境監(jiān)測、軍事、氣象、海洋、水文、地質(zhì)、冰川等很多方面有著廣闊的應(yīng)用[1]，并被認(rèn)為是最有效、最有潛力的衛(wèi)星傳感器之一。隨著多種衛(wèi)星SAR傳感器的成功發(fā)射，SAR應(yīng)用領(lǐng)域會更加廣闊。

但是，不同SAR傳感器產(chǎn)品的數(shù)據(jù)格式不同，造成了對SAR數(shù)據(jù)的應(yīng)用非常不便。GDAL（Geospatial Data Abstraction Library）[2]是一個在X/MIT許可協(xié)議下的開源柵格空間數(shù)據(jù)轉(zhuǎn)換庫，幾乎能夠讀取目前所有格式的遙感數(shù)據(jù)，很多著名的GIS類產(chǎn)品都使用了GDAL庫，包括ESRI的ArcGIS 9.2、Google Earth和跨平臺的GRASSGIS系統(tǒng)。因此在數(shù)據(jù)讀取方面GDAL是比較理想的選擇，但是它對屬性讀取的支持有限，所以單純用GDAL并不能滿足實際要求，還需要其他輔助庫。

1 數(shù)據(jù)的讀取和顯示

1.1 數(shù)據(jù)讀取

SAR數(shù)據(jù)讀取主要采用開源GDAL庫實現(xiàn)。默認(rèn)編譯的GDAL 能 夠 讀 取 ENVISAT、ERS、Radarsat、TerraSAR-X 的 數(shù)據(jù)格式產(chǎn)品，但是并不支持以HDF5格式保存的CosmoSkyMed產(chǎn)品，如果要支持HDF5數(shù)據(jù)格式，必須將HDF5庫編譯到GDAL中[3]。

在讀取SAR數(shù)據(jù)前，必須打開數(shù)據(jù)并獲取數(shù)據(jù)的數(shù)據(jù)集。不同的文件類型，數(shù)據(jù)集的獲取方法稍有差別。

對于非HDF5格式數(shù)據(jù)集的獲取方法：

GDALAllRegister（）； //注冊驅(qū)動

CPLSetConfigOption （“GDAL_FILENAME_IS_UTF8”，“NO”）； //支持中文路徑

GDALDataset*poDataset= （GDALDataset*）GDALOpen（lpszPathName， GA_ReadOnly）； //獲取數(shù)據(jù)集

poDataset即為獲取的數(shù)據(jù)集。

HDF5格式比較復(fù)雜，一個HDF5文件中可能包含一個或多個數(shù)據(jù)集。因此HDF5格式的數(shù)據(jù)集的獲取方法也相對比較復(fù)雜，下面是獲取所有數(shù)據(jù)集的方法[2]：

GDALRegister_HDF5（）；//注冊驅(qū)動

CPLSetConfigOption（“GDAL_FILENAME_IS_UTF8”，“NO”）； //支持中文路徑

GDALDataset*tempDataset；

vectordatasets； //數(shù) 據(jù) 集 類 型 的vector容器

GDALDataset*Dataset= （GDALDataset*）GDALOpen（lpszPathName，GA_ReadOnly）；

char** Ch_Sub_DS= Dataset->GetMetadata（“SUBDATASETS”）；

//獲取并保存HDF5文件中所有的數(shù)據(jù)集

if（ CSLCount（Ch_Sub_DS） >0 ）

{

for（int i=0；Ch_Sub_DS[i]!=NULL；i+=2 ）

{

string tmpstr=string（Ch_Sub_DS[i]）；

tmpstr=tmpstr.substr（tmpstr.find_first_of（“=”）+1）；

tempDataset=（GDALDataset*）GDALOpen（tmpstr.c_str（），GA_ReadOnly）；

datasets.push_back（tempDataset）；

}

}

datasets為 vector容器，里面存放著獲取的所有數(shù)據(jù)集。

由于1.8版本之后的GDAL默認(rèn)情況下不支持中文路徑，可以在編譯前修改，也可以在注冊完驅(qū)動后加上代碼：CPLSetConfigOption（“GDAL_FILENAME_IS_UTF8”，“NO”）[4]。

GDAL在SAR數(shù)據(jù)的讀取方面有著很大的優(yōu)勢，可以用統(tǒng)一的代碼讀取所有柵格數(shù)據(jù)，數(shù)據(jù)讀取的代碼有以下兩種方法：

第一種：GDALReadBlock （GDALRasterBandH Band， int xoff， int yoff， void*pData）。

第 二 種 ：CPLErr GDALDataset::RasterIO （GDALRWFlag eRWFlag， int nXOff， int nYOff， int nXSize， int nYSize， void* pData， int nBufXSize， int nBufYSize， GDALDataType eBufType， int nBandCount， int * panBandMap， int nPixelSpace， int nLineSpace， eBufType* nBufXSize， int nBandSpace）；

其中eRWFlag標(biāo)志是讀還是寫，如果是GF_Read則表示代碼為讀、如果是GF_Write則標(biāo)志代碼為寫；pData是存儲數(shù)據(jù)的按行存儲的一維數(shù)組。兩種方法各有優(yōu)缺點：第一種方法速度很快，但是每次讀取的行列數(shù)是固定的，由要讀取數(shù)據(jù)的存儲結(jié)構(gòu)決定，也不能抽樣讀取。第二種方法速度較慢，但是讀取方法靈活，不但能讀取任意的行列數(shù)，而且還能按照任意比例抽樣讀取。在實際編程中采用了第二種方法。因為絕大多數(shù)的圖像都是按行存儲，因此第二種方法按行讀取和存儲比按列要快的多。GDAL讀取數(shù)據(jù)上下顛倒，在圖像顯示時應(yīng)該注意。

1.2 數(shù)據(jù)顯示

圖像顯示用了另一個開源庫—OpenGL[5]。OpenGL顯示速度非?？於宜旧峡梢燥@示任意格式的像素，顯示的主要代碼：

void glDrawPixels （GLsizei width， GLsizei height，GLenum format， GLenum type， const GLvoid*pixel）；

其中format為像素格式，主要用到的有GL_LUMINANCE（灰度格式）和GL_RGB（RGB格式）；pixel是要顯示的數(shù)據(jù)的指針，也是按行存儲的一維或多維數(shù)組，與GDAL讀取的數(shù)據(jù)相對應(yīng)。另外合理使用OpenGL的雙緩存可以加快顯示速度，但并不是所有的顯卡都支持雙緩存，尤其是Intel的集成顯卡。

2 屬性的讀取

2.1 基本屬性的讀取

在使用SAR數(shù)據(jù)時，必須先獲取數(shù)據(jù)的一些基本信息，如，圖像的長、寬和投影信息等。在非HDF5的SAR數(shù)據(jù)格式中，GDAL可以用統(tǒng)一的代碼讀取不同數(shù)據(jù)格式的基本信息，如下所示：

GDALDataType GetRasterDataType（）； //獲取數(shù)據(jù)類型

int GetRasterXSize（void ）；//獲取圖像長

int GetRasterYSize（void ）；//獲取圖像寬

int GetRasterCount（void）；//獲取圖像一個數(shù)據(jù)集中波段數(shù)

GDALRasterBand*GetRasterBand（ int）； //獲取圖像中一個數(shù)據(jù)集中的某一個波段指針

virtual const char*GetProjectionRef（void）； //獲取圖像的投影信息

virtual CPLErr GetGeoTransform（double*）；//獲取圖像的仿射變換

virtual const GDAL_GCP*GetGCPs（）；//獲取圖像的控制點

對于HDF5格式的數(shù)據(jù)，GDAL支持不好，除長、寬、數(shù)據(jù)類型和波段數(shù)外，其它基本信息都無法正確獲取。因此必須借助其他的輔助庫，具體獲取方法在下節(jié)全部屬性的讀取中介紹。

2.2 全部屬性的讀取

GDAL對柵格數(shù)據(jù)屬性的支持并不好，即使是已經(jīng)編譯進HDF5庫的GDAL，也無法用統(tǒng)一的方法獲取其基本屬性，而且也不能讀取HDF5的全局屬性。非HDF5格式的數(shù)據(jù)也只能讀取很少一部分屬性。因此必須自己寫代碼或選擇其他的庫作為輔助。

2.2.1 ENVISAT和ERS全部屬性的讀取

ENVISAT和ERS不同SAR產(chǎn)品甚至相同SAR數(shù)據(jù)產(chǎn)品的不同模式，數(shù)據(jù)格式并不完全相同，如果自己編寫代碼將會比較麻煩，幸好屬性的讀取在其官網(wǎng)[6]都有源代碼，下載下來稍作修改就可以編譯成庫或者直接添加到軟件工程中使用。

2.2.2 Radarsat和TerraSAR-X全部屬性的讀取

Radarsat和TerraSAR-X都采用了XML格式，即數(shù)據(jù)放在Tiff數(shù)據(jù)中，屬性及Tiff路徑等重要信息放到XML中。兩種SAR數(shù)據(jù)所采用的XML與GDAL的XML域規(guī)范并不完全相同，為了能夠讀取XML中的屬性信息，選用開源的TinyXML庫，使用遞歸函數(shù)實現(xiàn)全部屬性的讀取。

2.2.3 CosmoSkyMed HDF5全部屬性的讀取

由于HDF5的自我描述性使其結(jié)構(gòu)異常復(fù)雜，導(dǎo)致GDAL對其支持不好，不僅全局屬性無法讀取，基本屬性也讀不全。為了能夠讀全其屬性，需要單獨編譯HDF5的庫，讀取屬性時用HDF5，讀取數(shù)據(jù)時用GDAL。需要注意的是GDAL也編譯進了HDF5的庫，而HDF5庫有些版本并不能向低版本兼容，因此，GDAL所用的HDF5庫和單獨編譯的HDF5庫版本最好相同。下面為以UTM投影為例，用HDF5庫讀取CosmoSkyMed HDF5格式投影信息的代碼：

char string_out[80]={0}；

hid_t attr， file；

hid_t atrr_type；

herr_t ret；

file=H5Fopen （lpszPathName， H5F_ACC_RDONLY，H5P_DEFAULT）； //打開 hdf5 文件

int nzone；

CString strProject；

//獲取投影方式

attr=H5Aopen_name（file，"Projection ID"）；

if（attr>=0）

{

atrr_type=H5Tcopy（H5T_C_S1）；

H5Tset_size（atrr_type， 50）；

ret=H5Aread（attr， atrr_type， string_out）；

strProject.Format（“%s”，string_out）；

ret=H5Tclose（atrr_type）；

ret=H5Aclose（attr）；

}

if（strProject== “UTM”）

{

//獲取投影帶號

attr=H5Aopen_name（file，"Map Projection Zone"）；

if（attr>=0）

{

ret=H5Aread（attr， H5T_NATIVE_INT， &nzone）；

ret=H5Aclose（attr）；

}

}

else//其他投影方式，在此略

{

…….

}

3 數(shù)據(jù)及屬性的保存

3.1 數(shù)據(jù)保存

數(shù)據(jù)保存的格式可以是GDAL支持的任意一種。與SAR數(shù)據(jù)的讀取一樣，保存前也需要獲取要保存的數(shù)據(jù)集，可通過下面3行代碼獲?。?/p>

const char*pszFormat= “GTiff”；

GDALDriver*poDriver=GetGDALDriverManager （）->GetDriverByName（pszFormat）；

GDALDataset *poDstDS = poDriver ->Create （pszDstFilename， width， height， n， eDT， NULL）；

其中pszFormat為要保存的格式；pszDstFilename為輸出路徑。為了后續(xù)操作的簡單，保存的格式應(yīng)該統(tǒng)一，考慮到格式的通用性，將數(shù)據(jù)保存成GeoTiff格式[7]。

對于數(shù)據(jù)的保存，GDAL也有統(tǒng)一的代碼，和讀取相同，保存也有兩種方法：第一種方法為GDALWriteBlock，參數(shù)和GDALReadBlock完全相同；第二中方法和數(shù)據(jù)讀取的第二種方法相同，只需要將參數(shù)eRWFlag改為GF_Write。

3.2 屬性保存

對于每種數(shù)據(jù)格式屬性數(shù)據(jù)都能保存到TXT中。但是由于有些屬性在后續(xù)處理中可能會用到，因此最好能保存以方便及時調(diào)用。GDAL可以方便的把屬性保存到XML域中，并容易讀出。

3.2.1 基本屬性保存

在基本屬性中，投影信息、仿射系數(shù)和控制點，不僅可以讀取也可以保存，下面是保存分別保存的代碼：

virtual CPLErr SetProjection （const char* ）；//保存投影信息

virtual CPLErr SetGeoTransform （double* ）；//保存仿射系數(shù)

virtual CPLErr SetGCPs （ int nGCPCount， const GDAL_GCP *pasGCPList， const char*pszGCPProjection ）； //控制點

接著2.2.3節(jié)中的例子，將從HDF5數(shù)據(jù)中讀出的投影信息保存，示例代碼：

OGRSpatialReference oSRS；

oSRS.SetWellKnownGeogCS（"WGS84"）；

oSRS.SetUTM（nzone）；

char*chproj=NULL；//投影信息

oSRS.exportToPrettyWkt（&chproj）；

poDataset->SetProjection（chproj）；

經(jīng)過這部操作后，只需要調(diào)用統(tǒng)一代碼virtual const char*GetProjectionRef（void），就可以獲取圖像的投影信息。仿射系數(shù)和控制點也可以采用類似的操作。

3.2.2 其余屬性保存

屬性的保存：

SetMetadataItem （const char*pszName，const char*pszValue， const char*pszDomain= “” ）；

其中pszName為屬性的名字，pszValue為屬性的值，pszDomain一般采用默認(rèn)值即可。

屬性保存后的讀取：

virtual const char*GetMetadataItem （const char*pszName， const char*pszDomain= “”）；

pszName為屬性的名字，和保存時采用的名字必須一樣，返回值即為屬性的值。

4 軟件實現(xiàn)

軟件采用MFC結(jié)合GDAL等開源庫實現(xiàn)了主要衛(wèi)星SAR數(shù)據(jù)的讀取、顯示、放大、縮小、1:1顯示、全圖顯示、區(qū)域選擇、漫游以及移動地理坐標(biāo)、圖像坐標(biāo)和像素值的顯示，以及圖像和屬性的保存。另外為了更方便的操作，軟件添加了鷹眼視圖、直方圖視圖和已打開文件列表視圖。以CosmoSkyMed hdf5格式數(shù)據(jù)為例，SAR圖像數(shù)據(jù)顯示如圖1所示，全部屬性保存如圖2所示。

圖1 SAR數(shù)據(jù)顯示Fig.1 SAR data display

圖2 SAR屬性保存示例Fig.2 SAR attributes save example

5 結(jié)束語

針對目前主要的衛(wèi)星SAR傳感器，文中研究了基于GDAL和HDF5對數(shù)據(jù)的讀取和保存；基于OpenGL數(shù)據(jù)顯示；以及基于GDAL、HDF5和TinyXML對屬性的讀取和保存，并著重介紹了GDAL和HDF5聯(lián)合讀取HDF5文件的方法。為了后續(xù)處理的簡便，最后將所有類型的數(shù)據(jù)轉(zhuǎn)化為Geotiff格式，并將屬性信息保存到與Geotiff文件名相同的XML文件中。軟件采用面向?qū)ο蠼Y(jié)構(gòu)，各功能用類實現(xiàn)，很容易在此基礎(chǔ)上添加其后續(xù)功能。對自主開發(fā)SAR處理軟件有一定的參考價值。

[1]Seelye Marti．海洋遙感導(dǎo)論[M]．蔣興偉，等，譯．北京：海洋出版社，2008：306-326．

[2]Norman Frank Warmerdam.GDAL（Geospatial DataAbstraction Library）.[EB/OL].[2012-04-20].http://www.gdal.org/.

[3]王繼成，蔣狄微，謝智劍．基于GDAL的HDF文件格式柵格數(shù)據(jù)提取的研究 [J]．計算機技術(shù)與信息發(fā)展，2011（8）：62-63．WANG Ji-cheng，JIANG Di-wei，XIE Zhi-jian.Research on GDAL-based HDF raster format data extraction [J].Computing Technology and Information Development，2011，（8）:62-63.

[4]李民錄．關(guān)于GDAL180中文路徑不能打開的問題分析與解決[EB/OL].（2011-07-16）[2012-05-12].http://blog.csdn.net/liminlu0314/article/details/6610069.

[5]Dave Shreiner， The Khronos OpenGL ARB Working Group．OpenGL編程指南[M]．7版．李軍，徐波等，譯．北京：機械工業(yè)出版社，2010．

[6]European Space Agency.ESA Observing the earth[EB/OL].[2012-05-22].http://www.esa.int/esaEO/index.html.

[7]牛芩濤，盛業(yè)華．GeoTIF圖像文件的數(shù)據(jù)存儲格式及讀寫[J]．四川測繪，2004，27（3）：105-108．NIU Qin-tao，SHENG Ye-hua.The storage and read/write of GeoTIFF image file[J].Surveying and Mapping of Sichuan，2004，27（3）:105-108.