衛星SAR數據讀取與保存方法研究與軟件實現

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

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

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

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

1 數據的讀取和顯示

1.1 數據讀取

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

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

對于非HDF5格式數據集的獲取方法：

GDALAllRegister（）； //注冊驅動

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

GDALDataset*poDataset= （GDALDataset*）GDALOpen（lpszPathName， GA_ReadOnly）； //獲取數據集

poDataset即為獲取的數據集。

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

GDALRegister_HDF5（）；//注冊驅動

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

GDALDataset*tempDataset；

vectordatasets； //數 據 集 類 型 的vector容器

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

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

//獲取并保存HDF5文件中所有的數據集

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容器，里面存放著獲取的所有數據集。

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

GDAL在SAR數據的讀取方面有著很大的優勢，可以用統一的代碼讀取所有柵格數據，數據讀取的代碼有以下兩種方法：

第一種：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標志是讀還是寫，如果是GF_Read則表示代碼為讀、如果是GF_Write則標志代碼為寫；pData是存儲數據的按行存儲的一維數組。兩種方法各有優缺點：第一種方法速度很快，但是每次讀取的行列數是固定的，由要讀取數據的存儲結構決定，也不能抽樣讀取。第二種方法速度較慢，但是讀取方法靈活，不但能讀取任意的行列數，而且還能按照任意比例抽樣讀取。在實際編程中采用了第二種方法。因為絕大多數的圖像都是按行存儲，因此第二種方法按行讀取和存儲比按列要快的多。GDAL讀取數據上下顛倒，在圖像顯示時應該注意。

1.2 數據顯示

圖像顯示用了另一個開源庫—OpenGL[5]。OpenGL顯示速度非常快而且它基本上可以顯示任意格式的像素，顯示的主要代碼：

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

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

2 屬性的讀取

2.1 基本屬性的讀取

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

GDALDataType GetRasterDataType（）； //獲取數據類型

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

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

int GetRasterCount（void）；//獲取圖像一個數據集中波段數

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

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

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

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

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

2.2 全部屬性的讀取

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

2.2.1 ENVISAT和ERS全部屬性的讀取

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

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

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

2.2.3 CosmoSkyMed HDF5全部屬性的讀取

由于HDF5的自我描述性使其結構異常復雜，導致GDAL對其支持不好，不僅全局屬性無法讀取，基本屬性也讀不全。為了能夠讀全其屬性，需要單獨編譯HDF5的庫，讀取屬性時用HDF5，讀取數據時用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 數據及屬性的保存

3.1 數據保存

數據保存的格式可以是GDAL支持的任意一種。與SAR數據的讀取一樣，保存前也需要獲取要保存的數據集，可通過下面3行代碼獲取：

const char*pszFormat= “GTiff”；

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

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

其中pszFormat為要保存的格式；pszDstFilename為輸出路徑。為了后續操作的簡單，保存的格式應該統一，考慮到格式的通用性，將數據保存成GeoTiff格式[7]。

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

3.2 屬性保存

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

3.2.1 基本屬性保存

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

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

virtual CPLErr SetGeoTransform （double* ）；//保存仿射系數

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

接著2.2.3節中的例子，將從HDF5數據中讀出的投影信息保存，示例代碼：

OGRSpatialReference oSRS；

oSRS.SetWellKnownGeogCS（"WGS84"）；

oSRS.SetUTM（nzone）；

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

oSRS.exportToPrettyWkt（&chproj）；

poDataset->SetProjection（chproj）；

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

3.2.2 其余屬性保存

屬性的保存：

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

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

屬性保存后的讀取：

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

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

4 軟件實現

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

圖1 SAR數據顯示Fig.1 SAR data display

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

5 結束語

針對目前主要的衛星SAR傳感器，文中研究了基于GDAL和HDF5對數據的讀取和保存；基于OpenGL數據顯示；以及基于GDAL、HDF5和TinyXML對屬性的讀取和保存，并著重介紹了GDAL和HDF5聯合讀取HDF5文件的方法。為了后續處理的簡便，最后將所有類型的數據轉化為Geotiff格式，并將屬性信息保存到與Geotiff文件名相同的XML文件中。軟件采用面向對象結構，各功能用類實現，很容易在此基礎上添加其后續功能。對自主開發SAR處理軟件有一定的參考價值。

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

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

[3]王繼成，蔣狄微，謝智劍．基于GDAL的HDF文件格式柵格數據提取的研究 [J]．計算機技術與信息發展，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]李民錄．關于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版．李軍，徐波等，譯．北京：機械工業出版社，2010．

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

[7]牛芩濤，盛業華．GeoTIF圖像文件的數據存儲格式及讀寫[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.