在实际应用中，空间数据种类繁多、关系复杂，且不同行业有不同内容的专业数据，空间数据的生产、维护分散在不同的单位且采用的数据网格标准各异，造成了部门或系统之间、不同历史阶段之间数据产品规格不统一^[1]。例如，卫星地面数据接收中心的原始遥感影像数据按照轨道景或轨道条带组织，而轨道景、轨道条带主要依据格网参考系GRS(Grid Reference System)和世界参考系WRS(Worldwide Reference System)所形成的固定参考网格并利用Path/Row对轨道景编码^{[2, 3]}；测绘部门的正射影像数据产品，高级数据产品(4级以上)按照地图图幅标准进行组织，低级数据产品按照轨道景数据组织；广泛使用的空间数据服务，如Worldwind、Google Earth、天地图等，均采用构建影像金字塔的方式进行数据组织。这种不同部门建立各自独立的遥感数据组织和索引方式，可很好地适应本部门内业务需要，但跨部门间数据应用时，不同的数据组织和记录方式给用户带来了很大的困难，不利于空间信息的检索、整合、共享与分发等。

为解决上述问题，本文提出基于GeoSOT (geographical coordinate subdividing grid with one dimension integer coding on 2ⁿ-tree)网格^[4]的全球遥感影像数据剖分组织，为遥感数据高效应用奠定基础。

近年来，国内外学者对全球离散网格^{[5, 6]}有很深入的研究，取得了很多成果，具体归纳为如下几类。第1类是基于柏拉图立体的正多面体(4面体、立方体、8面体、12面体和20面体)网格^{[7, 8]}，如O-QTM(octahedral-quaternary triangular mesh)^{[9, 10]}、SQT(sphere quadtree)^{[11, 12]}和STQIE(spherical triangle quadtree based on Icosahedron and ERLRP)^{[13, 14, 15]}；第2类是基于正六边形的全球离散网格^{[8, 16, 17, 18]}；第3类是基于Voronoi图的自适应全球离散网格^{[19, 20, 21]}；第4类是基于地理坐标划分的全球离散网格，如多级网格^[22]、Worldwind网格^[2]、Google Earth和Google Maps网格^[2]、Bing Maps网格^[23]、天地图网格^[23]和我国测绘数据网格^[24]等。

这些全球离散网格具有自身的特点。第1类全球离散网格可使得各网格单元在全球范围内的形状相似、各向同性，但各单元为球面三角形，从多面体到球体的映射关系比较复杂^[25]；第2类全球离散网格最紧凑、具有唯一的邻域，可用6个离散的速度向量描述连续的各向同性的流体^[26]，但上下层之间的网格单元不像三角形网格那样具有层次性和对称性，而是存在一定角度的旋转；第3类网格可针对空间数据的分布情况动态生成网格，适合组织空间点数据，但该网格单元不规则、无层次性，无法与比例尺和分辨率对应，且动态生成Voronoi图计算复杂，结构繁琐^[27]；第4类全球离散网格多是按经纬度网格进行划分，坐标转换方便，但高纬度地区存在一定的变形^[27]。针对遥感影像数据的矩形点阵结构的特点，基于三角形、六边形和Voronoi图等的网格，由于网格单元与影像数据块单元不一致，在组织管理遥感影像数据，需大量复杂数据转换工作，使得这3类网格仅用于理论研究和少数应用^{[28, 29]}，组织管理遥感影像数据均存在一定的局限性。基于经纬度坐标的全球离散网格单元与影像数据点阵相似以及坐标转换简单等优点，被Google Earth、Worldwind、天地图等广泛地应用。因此，本文主要对经纬度坐标全球离散网格讨论。同时，我国遥感影像数据产品多按地图图幅组织，Worldwind、Google Earth、Bing Maps和天地图等网格单元大小并不是整度、整分和整秒划分，不能无缝地拼接为地图图幅(具体见后面理论证明)，对跨部门之间遥感数据组织网格的兼容性和历史数据继承等问题考虑较少，在制作遥感影像数据产品还需数据拼接、裁剪等数据处理。

在此背景下，以北京大学程承旗教授为首的研究团队提出了GeoSOT剖分网格理论，期望利用GeoSOT网格独特的性质，解决跨部门之间、部门内各业务阶段遥感影像数据组织基准不统一的问题。

为了使得全球网格单元(剖分面片)大小为二进制整型、便于与地图图幅网格兼容，GeoSOT网格通过地球表面经纬度范围空间经过3次扩展后再对其进行严格的递归四叉剖分。具体GeoSOT网格划分方法如下：首先将地球通过简单投影变换到平面，将180°×360°地球表面空间扩展为512°×512°，并将该空间作为第0级剖分面片(cell)，且面片中心与赤道和本初子午线的交点重合(见图 1(a)和图 1(b))。在此基础上，对第0级剖分面片进行递归四叉剖分，直到1°剖分面片。为了进行整分剖分，对1°剖分面片先进行扩展，即从60′扩展为64′(见图 1(e))，再对扩展的剖分面片递归四叉剖分，直到1′大小的剖分面片。对于1′网格面片，同样先将其扩展到64″(见图 1(f))，然后递归四叉剖分直到32级(1/2048)″。这样，GeoSOT网格产生整度、整分、整秒和秒下级等多级剖分网格(见图 1(c))。这样，经纬度坐标空间经过3次扩展后，整个地球表面经纬度空间在经线方向和纬线方向可严格的整型二分，由此将整个地球分割为大到全球、小到厘米级的整度、整分和整秒的层次网格。

图 1 GeoSOT网格多级剖分与3次扩展示意图 Fig. 1 Multi-level subdivision and three extension of GeoSOT

图选项

针对GeoSOT网格剖分方案，采用64位编码对各级剖分面片进行标识。其中，第0级剖分面片标识为0，第1级剖分面片编码分别为00、01、02和03(见图 2(a))，并从第2级开始，下一级剖分面片的编码在上一级剖分面片编码后追加0、1、2和3。具体编码规则是，距赤道和本初子午线的交点最近的剖分面片为0，最远的为3，然后按照先沿纬线方向再沿经线方向对其他两个剖分面片分别为1和2(见图 2(b))。

图 2 GeoSOT网格单元的编码 Fig. 2 GeoSOT cells’ codes

图选项

由于GeoSOT网格是基于经纬度整度、整分、整秒的递归四叉剖分，与Google Earth、Worldwind、Bing Maps和天地图等格比较，GeoSOT网格与测绘数据地图图幅具有很好的同构性，证明如下。

为了方便后面的证明，首先介绍有关网格同构的几个概念，且讨论的全球经纬度网格是将球面网格投影到平面后的全球经纬度网格。

定义1：在全球经纬度网格G中，若映射f适合

点击浏览公式

式中，t、b、l、r代表网格单元的四边长度；θ代表任意两邻边的夹角，则称C为网格G的单元。其中，t、 b、 l、 r、 θ单位为度(°)。

第i层全球经纬度网格G的单元C，记作C(G，i)或C_i，表示如下

点击浏览公式

在全球经纬度网格G单元C_i=f(t_i，b_i，l_i，r_i，θ)中，参数取值如下

点击浏览公式

式中，L_t、L_b、B_l、B_r为整个地球或者第1层单元格的经度和纬度长度，单位为度(°)。当L_t、L_b、B_l、B_r、k₁、k₂、k₃、k₄和θ₀取不同值时，即可表示不同全球经纬度网格。

根据定义1，NASA Worldwind、Google Earth、Bing Maps、天地图、我国测绘数据规格网格和GeoSOT等网格单元C_i=f(t_i ,b_i ,l_i ,r_i ,θ)，其参数取值如表 1所示。

表 1 不同经纬度全球剖分网格单元f(t_i, b_i, l_i, r_i, θ)参数取值 Tab. 1 Values of f(t_i , b_i , l_i , r_i, θ) in different global subdivision grids

数据网格	ti	bi	li	ri	θ
Worldwind	ti=bi =360/(10×2i)，i∈N，Lt=Lb=360， k1= k2=10		li=ri =180/(5×2i)，i∈N，Bl=Br=180，k3= k4=5		90
Google Earth	ti=bi =360/2i，i∈N，Lt=Lb=360， k1= k2=1		li=ri =180/2i，i∈N，Bl=Br=180，k3= k4=1		90
Bing Maps	ti=bi =360/(2×2i)，i∈N，Lt=Lb=360， k1= k2=2		li=ri =180/(2×2i)，i∈N，Bl=Br=180，k3= k4=2		90
天地图	ti=bi =360/(2×2i)，i∈N，Lt=Lb=360， k1= k2=2		li=ri =180/2i，i∈N，Bl=Br=180，k3= k4=1		90
中国地图图幅					90
GeoSOT网格					90

表选项

定义2：对于全球经纬度网格G₁和G₂，C(G₁,i)= f₁(t_i，b_i，l_i，r_i，θ₁)和C(G₂，j)= f₂(t_j，b_j，l_j，r_j，θ₂)(i，j∈N)分别是G₁第i层和G₂第j层的网格单元。若n(n∈N)个相邻的C(G₁,i)覆盖的空间范围与C(G₂,j)覆盖的空间范围相同，则称G₁可聚合成G₂。特别的，G₁与G₂可以是同一个全球经纬度网格。

若t_j=λt_i，b_j=λb_i，l_j=μl_i，r_j=μr_i，θ₁=θ₂=θ且n=λμ(1≤λ≤n，1≤μ≤n，λ∈N，μ∈N，n∈N)时，称为长度聚合(如图 3)，表示为

点击浏览公式

或

点击浏览公式

图 3 长度聚合示意图 Fig. 3 Length of aggregation

图选项

定义3：对于全球经纬度网格G₁和G₂，C(G₁,i)(i∈N)和C(G₂ ,j)(j∈N)分别是G₁和G₂的单元。若有限个C(G₁,i)聚合而成的空间范围与C(G₂,j)的覆盖空间范围不相等，但形状相同，则称G₁弱聚合成G₂。

定义4：对于全球经纬度网格G₁、G₂和G₃，若G₁可聚合成G₂，而G₁可弱聚合成G₃，则称G₂与G₃弱同构。

因此，对于全球经纬度网格G₁和G₂，若G₁和G₂的θ相等，G₁和G₂必是弱同构。其证明见附录。

定义5：对于全球经纬度网格G₁、G₂和G₃，若G₁可聚合成G₂，同时G₁也可聚合成G₃，则称G₂与G₃同构。特别的，全球经纬度网格G与其本身同构。

根据定义4，若要证明GeoSOT网格与其他全球经纬度网格同构，只需证明GeoSOT网格可聚合成其他全球经纬度网格即可，具体证明如下。

证明：对GeoSOT网格G₁与其他全球经纬度网格G₂，假设C(G₁,i)(i∈N)和C(G₂,j)(j∈N)分别是G₁第i层和G₂第j层的网格单元。要证明G₁与G₂同构，也就是要证明存在正整数m(m∈N)，使得G₁中m个单元C(G₁,i)可聚合成G₂单元C(G₂,j)，即

点击浏览公式

根据定义2，证明式(6)成立等价于证明存在λ∈N，μ∈N，使得下面3个等式成立

点击浏览公式

点击浏览公式

点击浏览公式

而对式(7)，将表 1中的t_i代入得

点击浏览公式

点击浏览公式

因为k₁∈{1，15/16，225/256}，k₂∈{1，2，10，60，90}，因此式(11)变换为

点击浏览公式

在式(12)中，180/k₂∈N，k₁可以表示为k₁=15^α/2^2α，α∈{0，1，2}，并代入式(12)得

点击浏览公式

在式(13)中，180/k₂∈{2，3}，15^α≥1，只要使得2^i-j-8-β≥1即可，即

点击浏览公式

也就是，当GeoSOT网格层级i≥8时，总可使得式(13)中λ为正整数。

同理，将表 1中GeoSOT网格l_i和其他经纬度网格的l_j代入式(8)得

点击浏览公式

而k₃∈{1，15/16，225/256}，k₄∈{1，2，5，45，135/2}，则令k₃=15^α/2^2α，α∈{0，1，2}，k₄=2^β×3^γ×5^ε，β∈{-1，0，1}，γ∈{0，2，3}，ε∈{0，1}，则式(15)变换为

点击浏览公式

要使得μ为正整数，只要式(16)右边为正整数即可。而5^α+1-ε≥1，且当α和γ取特定的值时3^α+2-γ≥1，所以只要使得2^i-j-2α-β-7≥1，即

点击浏览公式

也就是，当GeoSOT网格G₁层级i≥6时，总可使得式(16)中μ为正整数。

因此，根据式(14)和式(17)，对于表 1中的其他全球经纬度网格G₂，只要GeoSOT网格G₁层级i≥7时，总存在正整数m=λμ使得式(6)成立，即证GeoSOT网格与其他全球经纬度网格同构。

而Worldwind、Google Earth、Bing Maps、天地图等网格与测绘数据规格网格具有弱同构性，可采用反证法进行证明。

证明：假设Worldwind、Google Earth、Bing Maps、天地图等网格与测绘数据规格网格同构性。

根据定义2，只要存在λ整数，使得式(18)成立，即证明Worldwind、Google Earth、Bing Maps、天地图等网格与测绘数据规格网格同构性

点击浏览公式

也就是

点击浏览公式

令j=i-q，则式(19)变换为

点击浏览公式

由于2^j始终为偶数，k₁={60，90}为2、3、5等素数的合数，因此要使得λ为整数，则式(20)中左边k必须为3的倍数。而Google Earth、Worldwind、Bing Maps、天地图等网格模型中k取值为{1,2,10}，不是3的倍数，所以式(20)中左边结果为非整数，即λ为小数，与假设λ为整数矛盾，所以假设为错误。而Worldwind、Google Earth等网格与测绘数据规格网格的夹角θ相等，因此Worldwind、Google Earth等网格与测绘数据规格网格弱同构性。

鉴于GeoSOT网格的多尺度性以及与我国测绘数据规格网格同构的特点，在不改变现有数据组织的前提下，本文提出通过逻辑剖分，建立以GeoSOT网格剖分面片为遥感影像景划分单元，实现基于空间区域的遥感数据“全球一张网”组织，提高多源多类型遥感数据使用效率。

我国经过几十年的发展积累了大量的遥感数据，若这些数据都按GeoSOT剖分面片进行物理分景处理，必是一个耗时耗力的巨大工程，且可能对目前遥感数据生产流程带来影响。因此，本文提出的基于GeoSOT遥感数据组织思路是，在不改变目前遥感数据组织模式的前提下，首先根据遥感数据的地理坐标范围，确定该数据的最小外包面片MBC(minimum bounding cell)，并将MBC面片编码作为整幅遥感数据的剖分标识，写入剖分索引表中，然后根据遥感数据的空间分辨率确定所在GeoSOT层，并通过GeoSOT剖分面片单元4个角点地理坐标与遥感数据的地理坐标对应，对遥感数据按GeoSOT网格剖分面片空间范围进行“虚拟划分”，即逻辑剖分，将遥感影像数据覆盖剖分面片的编码与剖分面片4个角点对应遥感影像数据像素行列号等信息，写入与遥感影像文件名相同的剖分索引文件中，由此建立基于GeoSOT剖分面片编码的遥感数据块剖分索引，从而实现遥感数据按空间区域数据组织。由于GeoSOT网格剖分面片覆盖全球且大小固定，与剖分面片关联的遥感数据也虚拟覆盖全球且遥感数据大小规格一定，形成遥感数据虚拟“全球一张网”的数据组织体系，如图 4所示。

图 4 遥感数据建立逻辑剖分的流程图 Fig. 4 Logical subdivision of remote sensing data

图选项

此方案中，遥感数据与GeoSOT网格层级l的对应关系采用的如下公式确定

点击浏览公式

式中，r为遥感数据分辨率，单位为°/像素；l为GeoSOT网格层级；k为常数，一般取2^-n，n∈N。

同时，由于各级剖分面片空间范围固定，而不同等级的影像数据产品几何精度不同，导致基于剖分面片“分割”的遥感影像数据会产生一些问题。例如，对地理精度误差较大的低级别(3级以下)的遥感影像，由于其存在较大的几何精度误差，导致本应该与剖分面片A关联的数据块，分割后可能与A相邻的剖分面片B关联。为了避免出现上述现象，在对不同地理精度误差的遥感影像数据逻辑剖分时，分别采用冗余逻辑剖分和精确逻辑剖分。对地理精度误差较大的低级别(3级以下)遥感影像数据逻辑剖分时，将剖分面片空间范围向周围扩展遥感影像地理误差的宽度，并按扩展后的剖分面片“分割”遥感影像，即为冗余逻辑剖分，从而使得地球表面与剖分面片地理范围对应区域的遥感影像数据始终在关联的剖分面片内。图 5为地理精度误差为d的遥感数据产品冗余逻辑剖分示意图，实线单元格为GeoSOT剖分面片，虚线网格为根据精度误差d扩展后的冗余剖分面片，冗余逻辑剖分即按虚线对影像数据进行“分割”。而对于地理精度误差较小的高级(4级以上)影像数据产品，直接根据GeoSOT剖分面片的空间范围进行逻辑剖分，建立遥感数据逻辑剖分索引，即为精确逻辑剖分，如图 6所示。

图 5 低级影像数据产品的冗余逻辑剖分示意图 Fig. 5 Redundant logical subdivision of lower-level remote sensing data products

图选项

图 6 高级影像数据产品的精确逻辑剖分示意图 Fig. 6 Precise logical subdivision of advanced remote sensing data products

图选项

基于GeoSOT遥感数据“全球一张网”的数据组织，使得遥感数据具有共同的数据组织基础，可避免不同部门、不同型号卫星遥感数据采用的数据组织网格不一致问题。同时使得大区域的遥感数据可直接聚合得到，即当利用遥感数据制作某规格数据产品时，根据规格数据产品的空间范围信息，例如地图图幅范围，计算出规格产品范围覆盖的剖分面片，并根据剖分面片编码提取相关的遥感数据。最后，根据GeoSOT网格与其他规格数据网格同构性，通过剖分面片聚合的方式即可快速生成规格数据产品(如图 7)，可有效减少参与计算的遥感数据数据量，从而提高遥感数据整合效率。

图 7 基于GeoSOT网格剖分面片聚合的数据产品生产过程 Fig. 7 Process of products aggregated by GeoSOT cells

图选项

根据GeoSOT网格及编码模型，在Worldwind开源代码的基础上，利用Microsoft Visual studio 2008 C#作为开发平台，Microsoft Office Access2007作为后端数据库系统，在CPU主频为2.4 GHz、内存为2 GB、硬盘转速为7200 rpm的PC机上进行了遥感数据GeoSOT组织与数据整合的试验验证。试验数据为局部地区某卫星1A、2B和3B级约32 GB遥感影像数据，分辨率分别为2 m、5 m和10 m，其格式为TIF文件。当式(21)中k取值为1/512，则上述分辨率的遥感影像分别对应GeoSOT网格第15层、第14层和第13层。对地理精度相对较低的1A级影像数据采用冗余逻辑剖分，而对地理精度相对较高的2B级和3B级影像数据采用精确逻辑剖分处理。3种不同级别影像建立逻辑剖分索引的耗时如表 2。

在影像逻辑剖分组织的基础上，利用GeoSOT网格与测绘地形图图幅之间的同构关系，通过剖分面片编码对单个剖分面片对应的影像数据以及不同区域、比例尺地图图幅的影像数据进行了数据检索和提取试验，其试验结果如表 3所示。

从表 3可看出，由于遥感数据建立了基于GeoSOT网格的逻辑剖分索引，在对剖分面片或规格产品范围对应的遥感数据进行数据整合具有较高的效率。

需要说明的是，本测试系统主要是验证遥感数据按GeoSOT网格进行逻辑剖分组织的效率和基于剖分面片聚合的数据整合的效率，以此说明基于GeoSOT的遥感数据组织的可行性和遥感数据应用的高效性。

本文对目前常用的全球经纬度网格和GeoSOT网格进行了分析，理论证明GeoSOT网格较Google Earth、Worldwind、BingMaps和天地图等网格，与现有的测绘数据规格网格具有更好的同构性。鉴于GeoSOT网格的多尺度性和对我国测绘数据规格网格同构性特点，在不改变现有遥感数据组织的前提下，提出了通过基于GeoSOT网格的逻辑剖分方法，实现基于空间区域的遥感数据“全球一张网”组织。同时，本文还提出了基于剖分面片聚合的规格数据产品生成方法，可有效减少参与计算的遥感数据量，从而提高遥感数据整合效率。最后，通过试验验证，通过对现有的遥感数据建立逻辑剖分索引，实施代价不大，且基于GeoSOT的统一数据组织有助于规格数据产品的数据整合。

全球经纬度网格G₁和G₂弱同构证明。

证明：对于全球经纬度网格G₁和G₂，C(G₁,i)= f₁(t_i，b_i,l_i,r_i，θ₁)和C(G₂，j)=f₂(t_j，b_j，l_j，r_j，θ₂)(i,j∈N)分别是G₁第i层和G₂第j层的网格单元，且t_j>t_i>0、l_j>l_i>0、θ₁=θ₂。

由于θ₁=θ₂，则C(G₁,i)和C(G₂，j)可以进行长度聚合。而要证明G₁和G₂弱同构，则只要证明存在正整数n₁和n₂，使得n₁t_i > t_j，n₂l_i>l_j即可。

取N₁=[t_j/t_i]+1，N₂=[l_j/l_i]+1，对于任意的正整数n₁和n₂，当n₁> N₁，n₂> N₂时，n₁t_i > N₁t_i > t_it_j/t_i= t_j，n₂l_i > N₂l_i >l_il_j/l_i= l_j同时成立。

即证，G₁与G₂是弱同构。