专利名称：一种面向虚拟地球的栅格数据快速建库方法
技术领域：
本发明涉及空间数据组织技术领域，尤其是涉及一种面向虚拟地球的栅格数据快速建库方法。
背景技术：
如今越来越多的空间数据以栅格的形式存储(包括地形数据和影像数据)。栅格数据成为了重要的数据源。随着对地观测技术的快速发展，人类能够快速地获取多时相、多波段、多分辨率的海量遥感影像数据。随着分辨率的不断提高，全球数据的数据量也在以惊人的速度增加，如仅福建省的30mX 30m分辨率的24位bmp格式的遥感影像就有700MB，若为 ImX Im分辨率，数据量将达到630GB。如何去存储、组织和管理栅格数据，如何高效地检索、 显示用户感兴趣的数据，使栅格数据能够更好地为实际应用服务越来越受人们的关注。有关文献高峰，安培浚.国际空间和对地观测技术发展战略新动向[J].遥感技术与应用，2008，(06) =686-696.；张更新，张昭，朱江.全球对地观测卫星现状及其产业发展综述[J].数字通信世界，2009，(10) 72-75.；王恩泉.中国版Google Earth的空间数据组织与管理研究[D].北京中国测绘科学研究院，2007.。随着Google Earth的发布与流行，虚拟地球技术进入了人们的视野。虚拟地球能够模拟显示大范围高分辨率的海量遥感影像数据，不仅用于各个行业领域内的专业应用研究，也能够满足人们足不出户即可感知整个地球。目前比较热门虚拟地球产品除了 Google Earth，还有美国航天航空局(NASA)的WorldWind、Bing 3D和Skyline等，国内不少公司与高校开展了相关研究，建立自主版权的虚拟数字地球平台，如武汉大学测绘遥感信息工程国家重点实验室自主研制的虚拟地球系统GeoGlobe。有关文献苗放，叶成名，刘瑞，孔祥生，李康荣，徐松浦，张远红.新一代数字地球平台与“数字中国”技术体系架构探讨[J].测绘科学，2007，(06) 157-158. ；http:// news. 3snews.net/industry/20081225/6054.shtml。虚拟地球通常采用多分辨率瓦片金字塔结构来组织栅格数据，以提供快速的网络调度和三维渲染。对栅格数据建库的过程也就是构建多分辨率金字塔的过程。通过构建多分辨率金字塔，在多尺度金字塔结构索引机制的配合下，可以快速获取指定分辨率和范围的栅格数据，并适合于网络传输和渲染处理，从而保证用户根据视景内的范围，在任一角度浏览的数据只是金字塔的子集。面向多分辨率瓦片金字塔的数据组织与调度是虚拟地球系统得以实现网络三维高效可视化的核心思想。如何对栅格数据快速构建多分辨率的金字塔一直是人们研究的热点。有关文献龚健雅，李斌等.当代GIS的若干理论与技术.武汉武汉测绘科技大学出版社.；吴晨琛.GeoGlobe中尺度空间数据集管理机制[D].武汉武汉大学，2008.。早期对金字塔的研究重点主要是如何提高成像质量和传输。1979年，Tanimoto等人提出了均值金字塔。均值金字塔的思想简单而直观，但应用到图像的渐进传输中却能获得很好的效果，引起了人们的关注。
有关文献朱程辉，何勇，王金玲.基于小波金字塔的快速图像匹配算法[J] 微计算机信息，2010，(10) 127-128.。同年，K. R. Sloan和Tanimoto又提出了 Reduced-sum金字塔的结构，该结构克服了均值金字塔数据量过大的问题，将金字塔中的像素数减少到与原始图像相同。有关文献刘国权，李守轩.基于小波图像金字塔的SSDA快速模板匹配算法[J]， 科技广场，2007，(11) 134-136.。但无论是均值金字塔还是Reduced-sum金字塔，其画面质量都不够理想，图像的混叠失真较大。于是在1983年，P. J. Burt等与J. Nao等提出了高斯-拉普拉斯金字塔结构，该结构在采样的同时利用高斯型低通滤波器对图像进行平滑，从而获得良好的视觉效^ ο有关文献丁素英.基于小波变换的Laplacian金字塔图像数据压缩[J].潍坊学院学报，2009，(04) =34-36.；叶勤，张小虎，王栋.一种基于区域生长与空间形状约束的高分辨率遥感图像道路提取方法[J].遥感信息，2010，(02) =25-29.。而后的学者对如何对大数据量的金字塔进行有效的存储和管理进行了研究。1989 年，又有学者提出了 Difference金字塔结构，通过将相邻层的像素相减来去除相关性，从而减少金字塔的信息熵。有关文献欧阳平，张玉方.形态学开闭运算在居民地边缘检测中的应用[J].测绘通报，2009，(01) :40-41·。随着计算机技术的快速发展，计算机能提供足够的空间存储金字塔数据，但是随着分辨率的提升，即使覆盖较小区域，一个遥感数据文件的尺寸也是比较大的，可能达到上十G，甚至上百G，从这样的超大数据文件(1个或者多个)构建出多分辨率瓦片金字塔，其瓶颈在于读操作。目前，一般采用的金字塔构建方法是逐层逐块构建金字塔，对每层金字塔都获取逐个像素的扫描行，然后将扫描行放人缓存，当缓存已满时采用邻近法替换数据。其缺点是多次扫描原始遥感文件，计算一个顶层瓦片时，涉及到原始遥感文件的较大范围，I/O效率低下。
发明内容
针对目前在虚拟地球系统中对栅格数据构建多分辨率金字塔的方法无法高效处理大数据量的遥感影像，为了提高构建多分辨率金字塔效率，本发明提供的技术方案是面向虚拟地球的一种针对海量栅格数据构建金字塔的方法，主要是针对tif、img等可以存储大数据量格式的文件。本发明的技术方案为一种面向虚拟地球的栅格数据快速建库方法，包括以下步骤
步骤1，读取原始遥感影像的信息，包括原始遥感影像的波段数bn、原始遥感影像的尺寸pwXph、原始遥感影像的影像分辨率pd以及原始遥感影像对应的地理范围(rw，rs, re, rn);定义瓦片结构体，所述瓦片结构体是用于在内存中存储瓦片数据并对瓦片进行标识的数据结构体，
步骤2，待构建金字塔的最顶层为分辨率最低的层次，标记为第0层；根据原始遥感影
5像对应的地理范围(rw，rs, re, rn)和第0层瓦片分辨率yd W]，计算出第0层瓦片的起止范围，计算方式如下，
sr
= (rw+180+le-8) /yd
er
= (re+180+le-8) /yd
sc
= (rs+90+le-8) /yd
ec
= (rn+90+le-8) /yd
其中，sr
为第0层的开始行、er
为第0层的结束行、sc
为第0层的开始列、 ec
为第0层的结束列、e为数学参数；
第0层瓦片分辨率yd
根据原始遥感影像的影像分辨率pd计算，计算公式如下， yd
=pdXtilesizeX2 厂1
其中，7为预设的待构建金字塔的层数，tilesize为预设的瓦片大小； 步骤3，比较第0层一个瓦片对应原始遥感影像所需的内存memvolumetO]与预设的可申请的内存，预设的可申请内存上限为ms个字节，当所需内存大小memvolumetO]不大于可申请内存上限ms时进入步骤4，否则进入步骤5 ；
所述第0层一个瓦片对应原始遥感影像所需的内存大小memvolumetO]根据原始遥感影像的波段数bn计算，计算公式如下， memvolume
=pixelvolume
Xbn
其中，pixelvolumetO]为第0层中存储一个瓦片覆盖原始遥感影像的其中一个波段的信息所需的内存大小，根据第0层瓦片分辨率ydW]、原始遥感影像的尺寸pwXph和原始遥感影像对应的地理范围(rw，rs, re, rn)计算，计算公式如下，
pixelvolume
= (yd
Xpw/(re-rw) +0. 5) X ( yd
X ph/(rn-rs)+0. 5) Xsizeof(type)
其中type为存储像素值的数据类型，sizeof ()是获取数据类型在计算机内占有字节数的函数；
步骤4，逐个处理第0层的每一个瓦片，处理方式为计算第0层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入到内存，记为数据块i， 数据块i对应的第0层瓦片及其四分后代瓦片的瓦片结构体记为集合Sti,将集合Sti初始化；遍历集合Sti,当遍历到集合Sti中的某个瓦片结构体时，在数据块i中获取相应的数据，通过插值获得瓦片数据并编码；
处理第0层的所有瓦片后进入步骤6 ；
步骤5，设某一层中存储一个瓦片对应的原始影像所需的内存大小为 1^1^01111^[1^61]，存储从第0层到该层上一层的所有瓦片结构体所需的内存大小是 tlstuvoIume [level],设 memvolume [level] +tlstuvoIume [level]取值最小的层为内存需求最小层wl ；进行以下步骤，
步骤5. 1，求取内存需求最小层wl，计算公式如下， wl=max (0, L(/-e_l)/2」)
其中e = Llog>/3)」，U代表向下取整沖表示第ο层瓦片个数，根据步骤2所得第0层的开始行sr
、第0层的结束行er
、第0层的开始列sc
和第0层的结束列ec
计算，计算公式如下，
6m=(er
-sr
+l) X (ec
-sc
+1)
步骤5. 2，将第0层到第wl-1层的所有瓦片结构体记为集合ts，将集合ts初始化，并加载到内存之中；
步骤5. 3，逐个处理第wl层的每一个瓦片，处理方式为计算第wl层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入内存，记为数据块 i，数据块i对应的第wl层瓦片及其四分后代瓦片的瓦片结构体记为集合Sti,将集合Sti初始化；遍历集合Sti和集合ts，
当遍历到集合Sti中的某个瓦片结构体时，在数据块i中获取相应的数据，通过插值获得瓦片数据并编码；
当遍历到集合ts中的某个瓦片结构体时，判断瓦片的覆盖范围是否包括数据块i，如果包括则插值获得该瓦片内相应部分区域的瓦片数据，在获得该瓦片内所有瓦片数据后进行编码；
处理第wl层的所有瓦片后进入步骤6 ； 步骤6，输出多分辨率金字塔。而且，步骤4和步骤5中，设置瓦片处理线程和瓦片计算线程两类线程，并开辟两块内存区间主缓存和备用缓存来存储数据块；
瓦片处理线程将原始遥感影像的数据块读入备用缓存，并在瓦片计算线程遍历完瓦片结构体后，将备用缓存设置为主缓存，并继续将下一个数据块读入备用缓存； 瓦片计算线程读取主缓存内的数据块，遍历瓦片结构体，插值计算瓦片数据。而且，瓦片结构体包括变量 count、status、level、row、col、pixelSpan 禾口指针 pTileBmp，其中变量count标识瓦片中有效值的个数、变量status标识瓦片状态、变量 level标识瓦片所在的层数、变量row标识瓦片的行号、变量col标识瓦片的列号、变量 PixelSpan标识瓦片中每个像素的分辨率，指针pTileBmp指向瓦片数据所在内存位置。本发明的优势在于从中间层(算法对于内存的需求是最小的层)开始构建金字塔， 打破了以往逐层构建的局限性，仅需扫描原始数据文件一次，减少了 1/0操作次数，消除了处理海量数据时读取数据的瓶颈问题，这对于超大文件来说提高了金字塔构建的效率。本发明还充分利用了现代计算机的多处理器多核，可并行计算的特点；采用瓦片处理线程和瓦片计算线程两类线程同时运行，充分利用了系统资源，通过配置多个瓦片计算线程，使程序读数据和写数据的速度相匹配，实现了读写操作的同步性，节省了读写操作相互等待的时间。
图1为本发明实施例的流程图2为本发明实施例的双线性内插示意图； 图3为本发明实施例的瓦片时空统一编码示意图。图4为本发明实施例的数据处理流程图。
具体实施例方式
为便于理解本发明技术方案起见，首先说明全球多分辨率金字塔全球多分辨率金字塔是基于球面的一种可以无限细分，但又不改变形状的地球体拟合格网，当细分到一定程度时，可以达到模拟地球表面的目的。本发明中采用的是通过等经纬度格网划分的全球多分辨率金字塔模型。在栅格数据金字塔中，按照全球等经纬度格网划分，每一层被裁切成统一形状、大小的数据块，这样的一个数据块称之为瓦片。这种全球多分辨率金字塔非常适合用来组织和管理多分辨率的海量卫星影像数据，因为它遵循以下规律
1)全球地理坐标的经度范围为[-180°，+180°]，纬度范围为[-90°，+90°]；
2)第k层金字塔的分辨率为第(k+Ι)层金字塔的2倍；
3)在任意层内，横向和纵向金字塔瓦片的数目之比为2 1 ；
4)金字塔瓦片的编码顺序由左到右，由下到上；
根据上述规则，很容易计算出每层瓦片的个数，同样也能快速计算出给定经纬度的点落在哪一个瓦片内。例如，已知一个点的经纬度为 (足劝，该点在第k层金字塔的行号RowNO和列号ColNo可以通过下面的公式计算出 RowNo= L ((0+ 90) / (180/2"")」mod f ColNo = L ((』+ 180) / (180/2"")」mod 2α+1) 其中U代表向下取整，mod是模运算。相反的，知道一个第k层的瓦片的行列号(X，Y)，同样也可以计算出瓦片的经度范围west east和讳度范围south north
west = ((χ mod 2(k+1)) X 180/2k) - 180 east = west + 180/2k south = ((y mod 2k) X 180/2k) - 90 north = south + 180/2k 因此，只需给出待构建金字塔的层次，每层的瓦片划分是已知的。本发明流程可采用计算机软件技术实现自动运行。以下结合附图和实施例详细说明本发明技术方案，参见图1，实施例流程如下
步骤1、首先进行预处理，包括统一空间参考、无效域去除和读元数据。实施例将WGS-84坐标系作为统一的空间基准，这样可以保证多源地理空间数据的无缝浏览，也提高了数据处理的效率。具体实施时，利用现有Erdas、Envi等遥感软件对遥感影像进行处理，即可将其投影坐标转换成WGS-84坐标系。无效域的去除有两种途径，即无效值过滤与多边形裁剪。无效值过滤方法比较简单高效，在需要去除成像不理想的区域时则可使用多边形裁剪。多边形生成可以采用手工勾勒方式，在手工勾勒的基础上还可以进行精细调位。实施例中，遥感影像成像效果较好， 仅仅影像边缘有一些无效区域，可以直接设定无效值，采用无效值过滤的方法进行无效域去除。实施例读取原始遥感影像(即元数据)的信息，通过调用Gdal开源库读取遥感影像的原始数据文件信息实现，包括原始遥感影像的波段数bn、原始遥感影像的尺寸pwXph、 原始遥感影像的影像分辨率Pd以及原始遥感影像对应的地理范围(rw，rs, re, rn)。实施例采用了现有技术中的 Gdal 库(Geospatial Data Abstraction Library)，Gdal 库是一个在X/MIT许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式，能够读取、写入、转换、处理各种栅格数据格式，并且是跨平台的。
实施例定义瓦片结构体 typedef struct tagMemTile{
int count;//定义count为整形变量 int status;//定义status为整形变量 int level;//定义level为整形变量 int row;//定义row为整形变量 int col;//定义col为整形变量 double pixelSpan;// 定义 pixelSpan 为浮数 FIBITMAP* pTileBmp;// 定义 pTiIeBmp 为指向 FIBITMAP 的指针 }MemTiIe;
其中tagMemTile为结构体名称，MemTile为结构体tagMemTile的变量名，count为瓦片中有效值的个数，level是瓦片所在的层数，row是瓦片的行号，col的瓦片的列号， pixelSpan是瓦片中一个像素的分辨率，FIBITMAP是FreeImage中用来存储瓦片数据的结构，pTileBmp是指向这个结构的指针；status为瓦片状态，可用于标识内存需求最小层的上层瓦片还是下层的四分瓦片以及瓦片是否已经处理完。步骤2，待构建金字塔的最顶层为分辨率最低的层次，标记为第O层；根据原始遥感影像对应的地理范围(rw，rs, re, rn)和第O层瓦片分辨率yd W]，计算出第O层瓦片的起止范围，计算方式如下，
sr
= (rw+180+le-8) /yd
er
= (re+180+le-8) /yd
sc
= (rs+90+le-8) /yd
ec
= (rn+90+le-8) /yd
其中，sr
为第O层的开始行、er
为第O层的结束行、sc
为第O层的开始列、 ec
为第O层的结束列、e为数学参数。实施例中，第O层瓦片分辨率yd
根据原始遥感影像的影像分辨率pd计算，计算公式如下，
yd
=pdXtilesizeX2 厂1
其中，7为预设的待构建金字塔的层数，tilesize为预设的瓦片大小，即一个瓦片包含的像素数目。步骤3，比较第O层一个瓦片对应原始遥感影像所需的内存memvolumetO]与预设的可申请的内存，预设的可申请内存上限为ms个字节，当所需内存大小memvolumetO]不大于可申请内存上限ms时进入步骤4，否则说明第O层中一个瓦片的相关数据不可以完全载入内存，进入步骤5。此时第O层瓦片的瓦片结构体的所需内存大小可以忽略。因此实施例中，第O层存储一个瓦片对应原始遥感影像所需的内存大小memvolumetO]根据原始遥感影像的波段数bn计算，计算公式如下，
memvolume
=pixelvolume
Xbn
其中，pixelvolumetO]为第0层中存储一个瓦片覆盖原始遥感影像的其中一个波段的信息所需的内存大小，根据第0层瓦片分辨率ydW]、原始遥感影像的尺寸pwXph和原始遥
9感影像对应的地理范围(rw，rs, re, rn)计算，计算公式如下，
pixelvolume
= (yd
Xpw/(re-rw)+0.5)X ( yd
Xph/(rn-rs)+0.5) Xsizeof(type)
其中type为存储像素值的数据类型，sizeof ()是获取数据类型在计算机内占有字节数的函数；
具体实施时，内存上限ms —般可设到500M，即500X1024X1024个字节。步骤4，逐个处理第0层的每一个瓦片，处理方式为计算第0层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入到内存，记为数据块i，数据块i对应的第0层瓦片及其四分后代瓦片的瓦片结构体记为集合Sti,将集合Sti 初始化；遍历集合Sti,当遍历到集合Sti中的某个瓦片结构体时，在数据块i中获取相应的数据，通过插值获得瓦片数据并编码；
处理第0层的所有瓦片后进入步骤6。为便于实施参考起见，提供实施例计算第0层中某瓦片在原始遥感影像的像素覆盖范围具体实现方法，设瓦片对应的地理范围为(m_RgnWest，m_RgnEast, m_RgnSouth, m_RgnN0rth)，该瓦片对应原始影像的像素范围为(xw，xe，ys，yn)，步骤1中读取遥感影像的元数据中遥感影像的左下角地面坐标为(x0，y0)，遥感影像的右上角地面坐标为 (xl，yl)，遥感影像像素在行方形上的地面分辨率为Dr，遥感影像像素在列方形上的地面分辨率为Dc，计算公式如下
Xw=( m_Rgnffest - x0 )/ Dc + 0. 5f ys= (yl - m_RgnSouth)/ Dr + 0. 5f xe= ( m_RgnEast - xO )/ Dc + 0. 5f yn=(yl - m_RgnNorth)/ Dr + 0. 5f 其中，f表示0.5是float类型。第0层为分辨率最低层，第0层一个瓦片对应遥感影像的数据量大于其他各层。因此如果第0层的一个瓦片可以读入内存，则只需逐个计算第0层的每一个瓦片及其所有四分瓦片，得到各瓦片中所有像素值，所得结果写入数据集中。步骤5，在第0层一个瓦片对应遥感影像不能完全载入内存时，实施例计算出所需内存容量最小的那一层，记为wl。即假设某一层中存储一个瓦片对应的原始遥感影像所需的内存大小为memvolumetlevel]，存储从第0层到该层上一层的所有瓦片结构体所需的内存大小是 tlstuvolume [level],设 memvoIume [level]+tlstuvoIume [level]取值最小的层为内存需求最小层wl。倘若实施时，发现memvolume [level]+tlstuvolume [level]仍大于内存上限ms，可以设计此时跳出程序，提示用户修改申请的缓存值。基于这种内存需求最小层设定，执行以下步骤
步骤5. 1，求取内存需求最小层wl。从第wl层开始计算多分辨率金字塔，算法对于内存的需求是最小的。这样可以根据设置的缓存大小确定从哪一层开始构建，实施例中Wl=Hiax(0，L a-e-l)/2」)
其中e = Llog>/3)」，U代表向下取整沖表示第ο层瓦片个数，根据步骤2所得第0层的开始行sr
、第0层的结束行er
、第0层的开始列sc
和第0层的结束列ec
计算，计算公式如下，m=(er
-sr
+l) X (ec
-sc
+1)
为便于理解本发明技术方案起见，下面给出了实施例中求取内存需求最小层层数方式的具体说明
设原始遥感影像有b个波段，输出的金字塔有/层，自上而下，第0层(分辨率最低层) 有《(!11彡1)个瓦片，第厂1层(分辨率最高层)则有个瓦片，总共有(m/3) X (47-1) 个瓦片。每个瓦片像素大小是&如果从第(/-1)层开始，即缓存第(/-1)层以上所有瓦片结构体，计为 fe/3) X(P^-I)，则所需内存为
Memvolume-I [7-1] +tlstuvolume U~l]=bXs2 X (Qn/3) X (4α_1)-1)) + bXs2X4° =bXs2(ipi/ ,) (4『"-Ι) + 4°)
如果从第(/-2)层开始，即缓存从第(/-2)层以上的所有瓦片，则所需内存为 memvolume [7-2] +tlstuvolume [7-2] =bXs2{ Qn/3) (4α_2)-1) + 41) 如果从第左(0彡左彡7-1)层开始，即缓存从第左层以上的所有瓦片，则所需内存为 Tk = memvolumeM+tlstuvolumeM=/ X52X (Qn/3) (Ai-I) + 4α+υ)。Tk+1_Tk = 3XbXs2X (Qn/3) Xa+2))，可见7；"和厶大小由Zli = / X(7+2) 来决定，其中，z = /3。设4e </K4(e+1)，由于 彡l，e彡-1，分两种情况来讨论
1) (he)彡2，此时，Zli彡4(时)_4时-2) ^=0，即‘恒为正值，因此左=0
时内存需求最小。2)a_e)>2，此时
设、l-e-2、=2h, Ah = nX4A-4(J-A_2) 彡 4(e+A)~4 (e+li=0, Ah^1 = / X< 4(e+A)-4(e+A+1)<0,即T她彡TkKT^iy因而々= 力时内存需求最小。设(/-e_2)=2A+l，zlA+1= /7χ4( +1)-4(Μ 彡 4(e+A+1)-4 (e+h)>Q,Ah = X 4A_4(7女2) < 4fe+A+1)-4(e+A+1)=0jpr_ 7；，因而仁0H1)时内存需求最小。综上所述，A=max(0， L (Pe-I)/2」)时，算法对内存的需求最小，其中e = Llog40 /3)」。当《=1，即第0层仅有1个瓦片，此时左=max(0，>^/2」)=1_//2」，即在等于>"//2」时内存需求最小。计算出内存需求最小层wl= max(0, L (/-e_l)/2」)后，就可以按对第0层到第 wl-1层、第wl层到第7 -1层分别按两种不同方式计算各层瓦片的瓦片数据。对第0层到第wl-1层，是将瓦片结构体初始化后加载到内存中，每从原始数据文件读出一个数据块， 遍历这些瓦片如果某瓦片中有像素所需数据位于该数据块内，插值计算出该像素的值。而对于第《1层及该层以下的瓦片则和步骤4处理方式一致。实施例采用步骤5. 2和步骤5. 3 实现。步骤5. 2，将第0层到第wl-1层的所有瓦片结构体记为集合ts，将集合ts初始化， 并加载到内存之中。例如第0层到第wl-1层某瓦片pMemTile的层行列为(k，rrl，ccl)，瓦片分辨率为Dc，瓦片大小为tilesize，则初始化该瓦片为
pMemTile=new MemTile ； //创建一个新的瓦片结构体 pMemTile_>count=0 ； //将 count 初始化为0 pMemTile->level=k;//对层号进行赋值pMemTile->row=rrl ；//对行号进行赋值 pMemTile->col=ccl;// 对列号进行赋值 pMemTile->pixelSpan=Dc/tilesize；// 计算瓦片的像素分辨率并赋值给pixelSpan pMemTi 1 e->pTi 1 eBmp=FreeImage_ Allocate (tilesize, tilesize, 32) ;//为 pTileBmp 分配内存 pMemTile_>status=l ；// 标识该瓦片为上层瓦片
步骤5. 3，逐个处理第wl层的每一个瓦片，处理方式为计算第wl层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入内存，记为数据块 i，数据块i对应的第wl层瓦片及其四分后代瓦片的瓦片结构体记为集合Sti,将集合Sti初始化；遍历集合Sti和集合ts，
当遍历到集合Sti中的瓦片结构体时，在数据块i中获取相应的数据，通过插值获得瓦片数据并编码；
当遍历到集合ts中的瓦片结构体时，判断瓦片的覆盖范围是否包括数据块i，如果包括则插值获得该瓦片内相应部分区域的瓦片数据，在获得该瓦片内所有瓦片数据后进行编码。当然，不包括则不进行任何处理。例如第wl层到第/ -1层某瓦片PMemTile的层行列为(k，rrl，ccl)，瓦片分辨率为Dc，瓦片大小为tilesize，则初始化该瓦片为
pMemTile=new MemTile ； //创建一个新的瓦片结构体 pMemTile_>count=0 ； //将 count 初始化为0 pMemTile->level=k;//对层号进行赋值pMemTile_>row=rrl ；//对行号进行赋值 pMemTile->col=ccl;//对列号进行赋值 pMemTile->pixelSpan=Dc/tilesize ；// 计算瓦片的像素分辨率并赋值给pixelSpan pMemTile_>pTiIeBmp=NULL ；//暂不为pTiIeBmp分配内存pMemTile->StatuS=-l;//标识该瓦片为wl层或下层瓦片
计算第层中某瓦片在原始遥感影像的像素覆盖范围，与步骤4计算第0层中某瓦片在原始遥感影像的像素覆盖范围具体实现方法类似
设第wl层某瓦片对应的地理范围为(m_RgnWest, m_RgnEast, m_RgnSouth, m_ RgnNorth),该瓦片对应原始影像的像素范围为(xw，xe, ys，yn)，步骤1中读取遥感影像的元数据中遥感影像的左下角地面坐标为(xO，yO)，遥感影像的右上角地面坐标为 (xl，yl)，遥感影像像素在行方形上的地面分辨率为Dr，遥感影像像素在列方形上的地面分辨率为Dc，计算公式如下
Xw=( m_Rgnffest - xO )/ Dc + 0. 5f ys= (yl - m_RgnSouth)/ Dr + 0. 5f xe=( m_RgnEast - xO )/ Dc + 0. 5f yn= (yl - m_RgnNorth)/ Dr + 0. 5f 其中，f表示0.5是float类型。当遍历完所有第wl层的瓦片后，即将原始影像扫描一遍，0到wl-Ι层的瓦片数据也就计算完毕(如果瓦片覆盖的范围超出了原始影像的范围，则以无效值处理)。对于集合 Sti中的瓦片结构体，通过插值获得瓦片数据后可以立即编码；对于集合ts中的瓦片结构体，因为每个数据块可能只对应瓦片内部分区域，可以将第wl层的所有瓦片对应的影像数据扫描一遍，即对原始文件扫描一遍后，再进行编码。为了便于提高瓦片的处理，具体实施时可以设定当0到wl-Ι层的某瓦片处理完后，将status的值设为0，之后便不用计算这个瓦片是否是在目前内存中数据块的范围内了。
为了提高处理效率，实施例在步骤4和步骤5中，设置瓦片处理线程和瓦片计算线程两类线程，并开辟两块内存区间主缓存和备用缓存来存储数据块；
瓦片处理线程将原始遥感影像的数据块读入备用缓存，并在瓦片计算线程遍历完瓦片结构体后，将备用缓存设置为主缓存，并继续将下一个数据块读入备用缓存； 瓦片计算线程读取主缓存内的数据块，遍历瓦片结构体，插值计算瓦片数据。根据硬件的条件不同，通过调整瓦片计算线程的数目，使程序读数据和写数据的速度相匹配，实现了读写操作的同步性，节省了读写操作相互等待的时间。步骤4和5中，具体插值获得瓦片数据实现可采用现有技术，为便于实施参考起见，提供实施例的插值处理方法说明如下
首先根据瓦片地理范围，计算其像素范围，将该范围的像素通过Gdal开源库读入内存。创建的金字塔瓦片里每个点不可能和全球基础数据文件中的每个像素一一对应，所以必须对其进行重采样。遥感影像采用最近邻采样算法，直接取与P(U)点位置最近像元F 的像素值为所需点的像素值作为采样值,即m=nm \为最邻近点，其影像坐标值为(、-,彻)
权利要求
1. 一种面向虚拟地球的栅格数据快速建库方法，其特征在于，包括以下步骤 步骤1，读取原始遥感影像的信息，包括原始遥感影像的波段数bn、原始遥感影像的尺寸pwXph、原始遥感影像的影像分辨率pd以及原始遥感影像对应的地理范围(rw，rs, re, rn)；定义瓦片结构体，所述瓦片结构体是用于在内存中存储瓦片数据并对瓦片进行标识的数据结构体；步骤2，待构建金字塔的最顶层为分辨率最低的层次，标记为第0层；根据原始遥感影像对应的地理范围(rw，rs, re, rn)和第0层瓦片分辨率yd W]，计算出第0层瓦片的起止范围，计算方式如下，sr
= (rw+180+le-8) /yd
er
= (re+180+le-8) /yd
sc
= (rs+90+le-8) /yd
ec
= (rn+90+le-8) /yd
其中，sr
为第0层的开始行、er
为第0层的结束行、sc
为第0层的开始列、 ec
为第0层的结束列、e为数学参数；第0层瓦片分辨率yd
根据原始遥感影像的影像分辨率pd计算，计算公式如下， yd
=pdXtilesizeX2 厂1其中，7为预设的待构建金字塔的层数，tilesize为预设的瓦片大小； 步骤3，比较第0层一个瓦片对应原始遥感影像所需的内存memvolumetO]与预设的可申请的内存，设预设的可申请内存上限为ms个字节，当所需内存大小memvolumetO]不大于可申请内存上限ms时进入步骤4，否则进入步骤5 ；所述第0层一个瓦片对应原始遥感影像所需的内存大小memvolumetO]根据原始遥感影像的波段数bn计算，计算公式如下， memvolume
=pixelvolume
Xbn其中，pixelvolumetO]为第0层中存储一个瓦片覆盖原始遥感影像的其中一个波段的信息所需的内存大小，根据第0层瓦片分辨率ydW]、原始遥感影像的尺寸pwXph和原始遥感影像对应的地理范围(rw，rs, re, rn)计算，计算公式如下，pixelvolume
= (yd
Xpw/(re-rw)+0.5)X ( yd
Xph/(rn-rs)+0.5) Xsizeof(type)其中type为存储像素值的数据类型，sizeof ()是获取数据类型在计算机内占有字节数的函数；步骤4，逐个处理第0层的每一个瓦片，处理方式为计算第0层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入到内存，记为数据块i， 数据块i对应的第0层瓦片及其四分后代瓦片的瓦片结构体记为集合Sti,将集合Sti初始化；遍历集合Sti,当遍历到集合Sti中的某个瓦片结构体时，在数据块i中获取相应的数据，通过插值获得瓦片数据并编码；处理第0层的所有瓦片后进入步骤6 ；步骤5，设某一层中存储一个瓦片对应的原始影像所需的内存大小为 1^1^01111^[1^61]，存储从第0层到该层上一层的所有瓦片结构体所需的内存大小是tlstuvoIume [level],设 memvoIume [level]+tlstuvoIume [level]取值最小的层为内存需求最小层wl ；进行以下步骤，步骤5. 1，求取内存需求最小层wl，计算公式如下，wl=max (O, L(/-e_l)/2」)其中e = Llog>/3)」，U代表向下取整沖表示第ο层瓦片个数，根据步骤2所得第O层的开始行sr
、第O层的结束行er
、第O层的开始列sc
和第O层的结束列ec
计算，计算公式如下，m= (er
-sr
+1) * (ec
-sc
+1)步骤5. 2，将第0层到第wl-1层的所有瓦片结构体记为集合ts，将集合ts初始化，并加载到内存之中；步骤5. 3，逐个处理第wl层的每一个瓦片，处理方式为计算第wl层中某瓦片在原始遥感影像的像素覆盖范围，并将该范围对应的原始影像数据的数据块读入内存，记为数据块 i，数据块i对应的第wl层瓦片及其四分后代瓦片的瓦片结构体记为集合Sti,将集合Sti初始化；遍历集合Sti和集合ts，当遍历到集合Sti中的某个瓦片结构体时，在数据块i中获取相应的数据，通过插值获得瓦片数据并编码；当遍历到集合ts中的某个瓦片结构体时，判断瓦片的覆盖范围是否包括数据块i，如果包括则插值获得该瓦片内相应部分区域的瓦片数据，在获得该瓦片内所有瓦片数据后进行编码；处理第wl层的所有瓦片后进入步骤6 ；步骤6，输出多分辨率金字塔。
2.如权利要求1所述面向虚拟地球的栅格数据快速建库方法，其特征在于步骤4和步骤5中，设置瓦片处理线程和瓦片计算线程两类线程，并开辟两块内存区间主缓存和备用缓存来存储数据块；瓦片处理线程将原始遥感影像的数据块读入备用缓存，并在瓦片计算线程遍历完瓦片结构体后，将备用缓存设置为主缓存，并继续将下一个数据块读入备用缓存；瓦片计算线程读取主缓存内的数据块，插值计算瓦片数据。
3.如权利要求1或2所述面向虚拟地球的栅格数据快速建库方法，其特征在于瓦片结构体包括变量 count、status、level、row、col、pixelSpan 禾口指针 pTileBmp，其中变量 count标识瓦片中有效值的个数、变量status标识瓦片状态、变量level标识瓦片所在的层数、变量row标识瓦片的行号、变量col标识瓦片的列号、变量pixelSpan标识瓦片中每个像素的分辨率，指针PTileBmp指向瓦片数据所在内存位置。
全文摘要
本发明涉及一种面向虚拟地球的栅格数据快速建库方法。该方法提出，如果第0层一个瓦片对应的原始遥感影像数据无法一次完全载入内存，从中间层(内存需求最小层)开始构建金字塔。这样打破了以往逐层构建的局限性，仅需扫描原始数据文件一次，减少了I/O操作次数，消除了处理海量数据时读取数据的瓶颈问题。本发明还设计了瓦片处理线程和瓦片计算线程两类线程，让两类线程同时运行，充分利用了系统资源，根据硬件不同的条件，配置多个瓦片计算线程，使程序读数据和写数据的速度相匹配，实现了读写操作的同步性，节省了读写操作相互等待的时间。提高了金字塔构建的效率，尤其适合于对超大遥感数据文件的处理。
文档编号G06F19/00GK102509022SQ20111036670
公开日2012年6月20日 申请日期2011年11月18日 优先权日2011年11月18日
发明者向隆刚, 曾珍 申请人:武汉大学