几种地球剖分网格分析与初探

当前的空间数据组织方面，尤其是地图数据，在统一组织方面已经形成一套完整的体系，如我国1992年颁布的《国家基本比例尺地形图分幅和编号GB/T130 89—92》标准等。针对遥感影像数据的组织与管理方面，在某类或者某系列卫星的影像数据统一组织上也有很完备的体系[1]，但缺少适合所有遥感卫星影像数据的统一组织上的网格框架，由此给基于遥感影像数据共享、检索、全球无缝表达、分发、编目、计算以及数据更新等应用带来诸多问题。本文从遥感影像数据特点出发，归纳了统一组织框架的设计原则和设计难点，分析了几种剖分格网作为遥感影像统一组织框架的优势和缺陷，为遥感影像统一组织提供参考。

遥感数据作为空间数据的一种，通常与传统的地图数据结合使用。因此，全球遥感影像的统一组织框架需考虑国家标准地图分幅体系与卫星数据分景体系之间兼容问题。

遥感卫星影像数据是关于地球表层或者地球表面附近区域信息的快照，具有与地球形状相似的性质。因此，全球遥感影像统一组织框架既要满足全球尺度无缝覆盖的要求，又要满足局部平面表达、精准量算的要求。如何在两者之间达到平衡，成为组织框架设计的难点之一。

我国现阶段对全球的空间数据采用经纬度进行组织，而对局部区域或者平面空间数据采用高斯坐标系下的地图分幅框架进行组织。鉴于我国空间数据组织现状，依据全球离散网格模型标准[2]，国内外研究比较多的地球剖分模型，如基于三角形或菱形的剖分模型[3-11]、基于正六边形剖分模型[12-16]、以及基于Voronoi图的剖分模型[17-20]，并不能很好地适合我国空间数据组织与管理。因此，在设计全球遥感影像数据统一组织框架时，在符合“Goodchild标准”[2]的前提下还需遵循如下原则：

(1) 可与地图分幅体系兼容。框架中各级不同大小的网格单元，又称剖分单元，能很方便地兼容或者成为地图分幅的有机部分，便于用于传统空间数据整合。

(2) 划分单元的形状尽量为矩形或正方形。框架中的网格单元形状最好为矩形或正方形，使其既可与影像中的像元点阵很好地对应，又便于影像的表达。

(3) 具有完整的编码体系。框架有一个完整的编码体系，使得大到地球、小到厘米级网格都有编码，且层与层之间的编码有很好的嵌套关系，以便于计算。

(4) 每个单元格具有明确的地理位置。框架中每个网格都有固定的空间范围，即每个网格中心、四个角点都有明确的坐标位置，可方便地与影像中坐标对应；每个网格编码在地球上具有准确的地理含义。

(5) 便于多源数据融合。在采用统一组织框架网格来切割影像时，可根据不同的分辨率，便于多源数据融合。

鉴于国内外学者对柏拉图立体、正多边形、Voronoi图等剖分模型进行了深入研究，并根据统一组织框架设计原则，本文只对等距、等纬差、正六面体、GeoSOT等几种剖分格网用于遥感影像统一组织的方案进行分析，具体如下。

等距剖分网格即基于CGCS2000大地坐标系统的地球椭球体，沿经线方向采用等距、沿纬线方向采用与纬圈相等的正方形建立全球等距格网框架。因此，等距剖分格网又分如下两种情况。

基于基点的等距剖分格网是指剖分单元以地球表面某点作为起始点形成全球无缝覆盖的等距剖分格网。比较典型的基点是地球两极点、本初子午线与赤道的交点、任意其他经线与赤道的交点。为了便于描述与计算，本文将其基点设为北极点，并假设正方形剖分单元的边长为d。基于北极的等距剖分格网示意图如图 1，图中剖分单元的边长与两个纬圈之间的距离相等，并采用不同的填充对不同距离的剖分面片进行区别。

基于CGCS2000椭球体[21-22]，若假设以正方形边长d为从北极点到某纬度之间沿大圆经线的弧长，则通过式 (1) 可算出该纬线圈所处的纬度，并由式 (2) 计算出纬圈的半径，由此得计算纬圈的周长。

式中，a为长半轴；e为第一偏心率；ϕ为纬度，E=(l0-d)/A0，l0为北极到赤道的距离；A0、R1、R2、R3见文献[23]中具体定义。

当d=25 km时，等距剖分格网分析统计如表 1所示。

由表 1可知，在第一圈剖分面片的周长与对应纬圈周长之差已经达到了42 917.863 m，并随着正方形单元从北极向赤道逼近，它们之间的差值越大。当d取不同的值时，也有类似的规律。因此，基于基点的等距剖分格网不能很好满足框架设计原则。

基于经线的等距剖分格网与基于基点的等距格网相似，只是在剖分面片的展开方式上，基于经线的等距剖分格网以某一经线作为“起跑线”，每一等距纬圈的剖分面片覆盖地球并由此展开 (如图 2)，剖分单元的边长与两个纬圈之间的距离相等，并采用不同的填充对不同纬圈的剖分面片进行区别。

采用式 (1) 可计算出不同纬圈的纬度，并由式 (2) 计算出各纬圈的半径。在此基础上，计算出各个纬圈内正方形面片的个数，即可计算在此高纬部分的各个正方形单元的重叠度，计算结果见表 2。

在每个纬圈内，正方形面片均不能完整地覆盖整个纬度带，且在高纬区域，越靠近极地，正方形单元之间的重叠度越大。因此，这种剖分方式也不能使得正方形无缝无叠地覆盖整个地球。

在传统的测绘数据管理上，多是采用等经差、等纬差方式对空间数据进行组织和管理，例如我国标准地形图分幅管理机制即采用这种方式。且当经差或者纬差小到一定程度时，其所代表地球表面的地理空间范围可以近似视为平面，如在赤道附近经度相差10′×10′的区域，其范围约为一个18.5 km×18.5 km。

若遥感影像统一组织框架采用纬差为φ′、以及纬差φ′对应的地表距离为正方形的边长为剖分单元进行等纬差剖分，其中，φ′=1、2、3、…、30。

基于式 (2) 的纬圈半径，可算出每一纬度所对应的纬圈周长、相邻两纬圈之间的距离，以及两者的比值，即为纬差φ′对应的正方形网格个数。为了计算方便，取φ′=10，其计算结果如表 3所示。

根据表 3结果分析，按此种剖分方案，某一地图图幅内剖分面片的分布如图 3。出现这种现象的原因是，由于从赤道到两极等纬差之间的地表距离并不相等，导致不同纬度带的正方形单元大小不相等，也不能完全部覆盖该纬度带对应的整个图幅。

针对剖分单元为矩形或正方形的原则，在柏拉图立体中，正六面体内接地球，有可能使得地球表面剖分为由多个四边形组成，如图 4。

基于图 4中正六面体的剖分框架具体剖分步骤如下：

(1) 设地球球面内接的正六面体棱长为d，赤道面与六面体的一面平行。六面体对角线等于球面直径，即，因此六面体的8个顶点经纬度分别为

即8点为 (35.264 390°，0°)、(-35.264 390°，0°)、(35.264 390°，90°)、(-35.264 390°，90°)、(35.264 390°，180°)、(-35.264 390°，180°)、(35.264 390°，-90°)、(-35.264 390°，-90°)，其相应的直角坐标容易计算，此处略。

(2) 从步骤 (1) 任选一面进行讨论，不妨选如下四边形：(35.264 390°，0°)、(35.264 390°，90°)、(-35.264 390°，90°)、(-35.264 390°，0°)。由此可知0级剖分对应的边长为赤道周长的四分之一，即约10 008 014.6 m。

(3) 将上一步得到的四边形，逐个进行 (近似)4等分：分别取各边中点，并用测地线 (大圆弧) 连接对边的中点 (如图 4(c))。

(4) 对同一分辨率的每个四边形，重复步骤 (3)。

(5) 当分辨率逐级提高 (四边形边长逐级减半) 时，对上一级的四边形逐个重复步骤 (3)、(4)。

由此即可得到由不同层级、不同尺度的剖分面片无缝覆盖地球表面的剖分体系。

下面对基于正六面体剖分框架中的剖分面片的面积变形、边长变形和角度变形进行分析。

该剖分体系中的剖分面片在地球表面为一个椭球面四边形 (如图 5)。设任意一个椭球面四边形ABCD，其面积的计算方法：

先用余弦公式求出四边形的4个球面角，然后代入式 (3) 即可得到面积：

式中，A、B、C、D分别为椭球四边形ABCD顶点所对应的内侧球面角。球面角的计算方法：

设AD、AB分别为过A点沿着曲线AD、AB的切线，相应的切面法向量分别为nAD、nAB，则球面角∠DAB(即角A) 满足

通过式 (3) 和式 (4)，即可分析出不同剖分层级上剖分面片的面积变形情况，如图 6和表 4所示。

由图 6可知，以初始剖分 (0级) 中心线为参照 (下同)，沿水平方向，剖分面片的面积变形随着经度的增加而增加；沿竖直方向，随着纬度的增加而增加。可归纳为：距离初始剖分中心线越远，剖分面片的面积变形越大。

剖分面片的边长，等于过该球面曲线两端点大圆的劣弧之长，用球面距离公式即可得到。图 7和表 5是部分等级剖分面片的边长变形分析结果。

剖分面片内侧夹角，即为该四边形顶点所对应的球面角，计算方法见式 (4)。图 8和表 6是部分等级剖分面片的角度变形分析结果。

根据上述分析可知，采用正六面体剖分地球的方式，虽然有可能到一定尺度后的剖分面片为正方形、且不同层级之间的剖分面片具有较好的递归关系，但存在如下问题：

(1) 虽然正六面体的6个面为正方形，但根据其等分递归生成的剖分面片并不都是正方形，只是特定的层级上的剖分面片接近正方形；

(2) 剖分面片与传统测绘中的经纬度并没有很明确的关系，尤其是两极剖分面片，如图 4(a)中的⑤、⑥对应的区域。因为在0级剖分面片⑤和⑥区域，在细分到1级剖分面片时，虽然可以与地球±45°和±135°经线重合，然而在继续按四等分剖分时，不再可能与经线、纬线重合。因此会增加与传统测绘数据之间转换的难度。

(3) 由于正六面体的8个顶点所在地球表面的纬度比较特殊，南北纬35.264 390°，则在与地图分幅的对应关系上，必然会割裂该区域的地图图幅。

因此，基于正六面体的格网剖分框架并不能很好地满足遥感数据统一组织框架设计原则。

为了解决跨部门间、部门内各业务阶段遥感影像数据组织基准不统一的问题，以北京大学程承旗教授为首的研究团队提出了GeoSOT剖分网格[24-26]。GeoSOT网格首先将地球表面空间扩展为512°×512°，面片中心与赤道和本初子午线的交点重合。在此基础上递归四叉剖分至1°剖分网格。然后将1°网格扩展为64′，再对其递归四叉剖分至1′剖分网格。最后将1′剖分网格扩展至64″，再递归四叉剖分至 (1/2048)″(如图 9)。在此基础上利用64 bit的反Z序对称编码对每个剖分面片进行唯一标识[24](如图 10)。

通过GeoSOT划分框架可知，GeoSOT网格面片不仅可实现全球覆盖而且可对每个剖分面片进行唯一标识。但由于GeoSOT网格是基于等经纬度进行划分，GeoSOT网格剖分面片没有角度变形，其边长长度沿经线方向没有变形 (图 11(b))、沿纬向方向随着纬度增加逐渐变小 (图 11(a))，剖分面片的面积也是随着纬度的增加逐渐变小 (图 12)。其中，剖分面片面积计算公式与式 (3) 相同。

且经过研究证明，GeoSOT网格与Google Earth、Worldwind、Bing Maps和天地图等格比较，GeoSOT网格与测绘数据地图图幅具有很好的同构性[24]，更适合对传统空间数据和遥感数据进行组织管理。

综上，GeoSOT网格符合全球遥感数据统一组织框架设计的各项原则，可用于遥感影像的统一组织。

通过对等距、等纬差、正六面体等几种剖分格网方案分析，采用矩形或正方形理论上不能无缝覆盖地球表面或者其剖分面片编码与传统测绘数据之间进行快速地转换，而GeoSOT剖分网格能够很好地与我国测绘数据地图图幅框架同构，可作为全球遥感影像统一组织框架。