全球空间网格将地球空间离散化为各个剖分等级的网格单元，为每一个网格单元建立全球唯一的网格编码，以网格单元为基础实现全球空间信息的统一组织，为全球范围的空间数据管理、空间信息分析、空间问题建模提供统一框架^[1]。按照网格剖分对象，全球空间网格可分为球面网格和球体网格。常见的球面网格包括球面正多面体网格^[2]、球面Voronoi网格^[3]、混合式球面网格^[4]等。常见的球体网格包括阴阳网格^[5]、立方体网格^[6]、球体退化八叉树网格^[7]、圈层立体网格^[8]等。

利用空间填充曲线的降维和聚簇性，设计高效的网格单元编码算法，是全球离散网格领域的一个重要研究方向^[1]。文献[9]提出了基于Hilbert曲线的全球离散网格影像金字塔索引算法，文献[10]基于Z曲线设计了球体退化八叉树网格的编码算法，文献[8]提出基于Hilbert曲线的扩展八叉树剖分网格编码算法，文献[11]将Hilbert曲线扩展应用到地月一体圈层立体网格的编码算法中。Hilbert曲线编码的经典方法是由文献[12]针对二维空间构造提出的。该算法是基于Morton码的二进制位操作法，算法复杂度为O(n²)。文献[13]根据Hilbert曲线的空间层次分解特征，通过栅格空间层次分解与构造域状态转移向量，设计了二维空间中Hilbert码的递归生成算法，算法复杂度为O(max (n))。文献[14]设计了基于状态转移矩阵的Hilbert码快速生成算法，将Hilbert状态转移矩阵转换为C++中的数组运算, 减少了Hilbert码计算过程中的嵌套循环及迭代处理，将算法复杂度降为O(n)。上述编码方法都只是在二维空间中进行讨论的。文献[15]研究了n维Hilbert曲线的生成问题，提出基于静态演化规则、自底向上的生成算法，算法复杂度为O(nlog n)。文献[16]设计了任意维空间中具有线性复杂度的Hilbert码算法。上述Hilbert曲线研究工作主要集中在曲线生成与排列码计算方面，对全球离散网格数据索引设计与优化发挥重要作用。但是，以Hilbert码为基础的网格编码计算分析研究却鲜见报导。

随着网格剖分模型研究的深入，全球空间网格已经在地球系统模拟器、全球气候变化、公众地理信息服务等科学研究和经济社会的诸多领域得到应用。但是，全球空间网格不仅仅是一种空间数据组织管理框架，而且将发展成为新一代空间分析模型。网格编码代数是这一发展过程中重要的理论推动^[17]。在网格编码上构建代数运算法则，将传统上复杂耗时的空间坐标计算转换为网格编码的代数运算，不仅给很多算法带来效率上的显著提升，还将为空间分析带来新的模型和方法。这一点已经在遥感图像处理^[17]等领域出露端倪。

本文首先基于Hilbert曲线设计八叉树立体网格单元编码，然后研究Hilbert曲线层级演进关系的计算方法，最后给出用于网格编码代数运算的Hilbert曲线网格编码若干操作算子。

1 八叉树立体网格Hilbert编码

对于八叉树结构的立体网格，可采用“行向编号-列向编号-径向编号”结构设计网格单元编码，简称“行列高”(I, J, K) 码。若剖分层次为N，则八叉树立体网格的大小为2^N×2^N×2^N，可对其中的每一个网格单元赋予3维N阶Hilbert编码作为唯一网格编码。图 1给出了3维1阶、2阶Hilbert编码与网格空间的对应关系。表 1所示为3维1阶Hilbert编码顺序及其对应二进制码与空间坐标 (I, J, K) 关系。

3维2阶Hilbert曲线是将3维1阶Hilbert曲线的节点网格空间进行再次细分得到的，如图 1(b)所示，其对应得编码等于“1阶编码+细分编码”。例如3维1阶Hilbert曲线序列码为3(011) 的网格被细分后获得对应3维2阶Hilbert编码如表 2所示。

根据3维1阶、2阶Hilbert编码生成过程，可以从图形结构上发现Hilbert编码的两个层级之间的关系为

式中，H_N+1²表示3维N+1阶Hilbert编码的二进制形式；H_N²表示3维N阶Hilbert编码的二进制形式；H₁^2′表示该立体网格单元细分1次后的子单元Hilbert编码 (二进制形式)。因此，Hilbert曲线序列码随剖分层级的变化中包含了父级与子级网格单元之间的空间关系，为基于网格编码构建代数运算法则，进而实现网格空间分析提供了可能的途径。

2 Hilbert曲线的层级演进关系

下面在分析Hilbert曲线生成过程中，由Hilbert曲线数学性质推导分析N阶与N+1阶Hilbert曲线序列码之间的关系。

2.1 Hilbert曲线地址

编码空间限定为M(M≥2) 维欧拉空间R^M，与其对应的是M维N阶Hilbert曲线H_N^M。

用S_a^{N, M}表示一个M维超立方体，其中N为剖分等级，a表示该超立方体的地址。初始超立方体未进行剖分，则N=0，a=0，记为S₀^{0, M}。在每个维度上对S₀^{0, M}进行N次二等分，得到子超立方体集合{S_a^{N, M}}。设S₀^{0, M}各个边长相等且均为单位长度1，则S₀^{0, M}和S_a^{N, M}的体积分别为1和 (1/2^N)^M，即是S₀^{0, M}由2^NM个S_a^{N, M}组成。

考虑到M维欧拉空间R^M与1维Hilbert曲线之间映射是一一映射，则有

Hilbert曲线具有自相似的层次结构，因此当把某一个S_a^{N, M}再细分1次时，得到2^M个子单元，每个子单元的地址可由两部分二进制段衔接起来，第1部分是a，第2部分是第i个子单元在S_a^{N, M}中的相对位置，记为b_i。可得

式中，b₁, b₂, …, b_2^M是一个M位二进制数，且满足当i≠j时b_i≠b_j。

2.2 Hilbert曲线起点和终点

Hilbert曲线起点和终点必须经过“垂直”子超立方体。“垂直”子超立方体是指其中包含S_a^{0, M}的任意一个顶点的子超立方体。集合{S_a^{N, M}}中有且仅有2^M个“垂直”子超立方体，每个“垂直”子超立方体都可以作为Hilbert曲线的起始点。对于每一个起始点，存在M个对应的终点，且起点和终点是相互邻近的。因此，Hilbert曲线所有起始点和终点组合的数量是M×2^M。

“垂直”子超立方体可被分为A和B两类。A类是指包含M个0的一组，其余为B类。由于Hilbert曲线的起点和终点是相互邻近的，若起点在A类中，则终点在B类中，反之亦然。例如，当M=3时，一条Hilbert曲线的起点和终点可以是A={000, 011, 110, 101}，B={001, 010, 111, 100}。因此，在研究Hilbert编码层级关系过程中，只关心起点位于A类中的Hilbert曲线，故Hilbert曲线数量减少至M×2^M-1。

2.3 Hilbert曲线地址列表

为了标识不同起点终点的Hilbert曲线，使用φ_l^{N, M}(h) 表示第l条Hilbert曲线序列码h对应的S_a^{N, M}的地址，其中l=1, 2, …, M×2^M-1。

M维空间第1条1阶Hilbert曲线上各个地址形成的列表φ₁^{1, M}可以从式 (4) 所示的递推公式计算得到

式中，0φ₁^{1, M-1}表示二进制“0”与φ₁^{1, M-1}的连接，φ₁^{1, M-1}表示对应的二进制位反转 (即0变为1, 1变为0)。例如，φ₁^{1, 3}={000, 001, 011, 010, 110, 111, 101, 100}是3维1阶第1条Hilbert曲线对应的地址列表。

下面将分3步确定其他Hilbert曲线对应的地址码φ_l^{1, M}(l=2, …, M2^M-1) 列表。

(1) 确定起始点地址

Hilbert曲线集合φ_l^{N, M}(h)(l=1, 2，…, M2^M-1) 符合性质1；φ_l^{1, M}(1) 和φ_l^{1, M}(2^M) 的汉明距离 (二进制位不同的个数) 等于1，且φ_l^{1, M}(1) 和φ_l′^{1, M}(1)(l≠l′) 的汉明距离必然是偶数^[18]。

为了确定哪些φ_l^{1, M}(l=2, …, M2^M-1) 满足上述性质，将Hilbert曲线序列码H为奇数的φ_l^{N, M}(h) 地址列表组成一个集合，例如φ₁^{1, M}(1), φ₁^{1, M}(3), …, φ₁^{1, M}(2^M-1)。由于它们之间的汉明距离是偶数，因此将其放到一起构成一个集合，即前面所说的类型A。考虑到一个Hilbert曲线起点对应的终点数量是M，可以使φ₂^{1, M}(1), φ₃^{1, M}(1), …, φ_M^{1, M}(1) 表示与φ₁^{1, M}(1) 具有相同起点。同样，φ_M(s-1)+1^{1, M}(1), φ_M(s-1)+2^{1, M}(1), …, φ_M(s-1)+M^{1, M}(1) 具有与φ₁^{1, M}(2s-1) 相同的起点，其中s=2, 3, …, 2^M-1。这样就能确定所有Hilbert曲线的起始点地址。

(2) 确定终点地址

φ21, M(2M), …, φM1, M(2M), …, φM2M-11, M(2M) 的终点地址可以从各自的起点地址计算得到。对于每一个起始点地址φ11, M(1)=φ21, M(1)=…=φM1, M(1)，其与相应的终点地址汉明距离等于1。例如，φ11, M(1)=000…0(M位二进制数，都是0)，则对应的终点地址为100…0，010…0，001…0，…，000…1(M位二进制数，只有一位是1，其他均是0)。将这些终点地址分别指定给φ11, M(2M), φ21, M(2M), …, φM1, M(2M)。总之，可以通过对M位二进制数的一个二进制位进行置反计算对应的终点地址。定义该运算函数为InvertBit (a, m)，其含义是将a对应的地址 (M位二进制数) 的第m位 (从高位到低位计数) 置反。因此，可将终点地址计算用如下函数表达

式中，s=1, 2, 3, …, 2^M-1；m=1, 2, …, M。

(3) 确定起点和终点之间其他地址的顺序

起点地址集合{φ_l^{1, M}(1)}和终点地址集合{φ_l^{1, M}(2^M)}中的对应元素分别作为第l条Hilbert曲线的起点和终点，则第l条Hilbert曲线剩余点的地址集合为{φ_l^{1, M}(h)|h=2, 3, …, 2^M-1}。

下面以2维和3维网格空间为例 (图 2)，讨论给定Hilbert曲线起点和终点情况下剩余点的地址顺序集合。在2维网格空间中，若指定00和01分别是Hilbert曲线的起点和终点，则其对应的Hilbert曲线路径只有一条 (00, 10, 11, 01)。在3维网格空间中，若指定000和100分别是Hilbert曲线的起点和终点，则其对应的Hilbert曲线路径有4条，分别为 (000, 001, 011, 010, 110, 111, 101, 100)、(000, 001, 101, 111, 011, 010, 110, 100)、(000, 010, 011, 001, 101, 111, 110, 100) 和 (000, 010, 110, 111, 011, 001, 101, 100)。因此，针对每一组起点终点组合，需要选定一组地址顺序列表，并将其指定给{φ_l^{1, M}(h)}。

基于φ₁^{1, M}的地址顺序列表，可采用交换X₁和X_l(l=2, 3, …, M) 两个坐标轴的方法获得的φ₁^{1, M}地址顺序列表。若定义位操作函数ExtractBit (φ, M) 是指提取二进制数a的第m位，则φ_l^{1, M}地址顺序列表的生成操作可以表示为交换ExtractBit (φ₁^{1, M}(h), 1) 和ExtractBit (φ₁^{1, M}(h), l)。因此，生成φ_l^{1, M}(l=2, 3, …, M) 地址顺序列表的公式为

式中，h=2, 3, …, 2^M-1是指Hilbert曲线序列码；φ₁^{1, M}(2^M) 和φ_l^{1, M}(2^M) 是第1条和第l条Hilbert曲线的终点地址；⊕表示二进制异或运算符。

同样地，需要确定{φ_M+i^{1, M}}(i=1, 2, 3, …, M) 的地址顺序列表，其可以根据Hilbert曲线如下的性质2，从{φ_i^{1, M}}计算得到。

Hilbert曲线符合性质2^[18]：给定一个地址顺序列表φ_i^{1, M}和一个起点地址φ_i′^{1, M}(1)(φ_i′^{1, M}(1)≠φ_i^{1, M}(1))，则新的地址顺序列表φ_i′^{1, M}将从φ_i′^{1, M}(h+1)=φ_i′^{1, M}(h)⊕(φ_i^{1, M}(h)⊕φ_i^{1, M}(h+1))(h=1, 2, …, 2^M-1) 依次计算得到^[18]。

以φ_i^{1, M}为基础，可以计算φ_M+i^{1, M}的地址顺序列表。重复该过程，可以根据φ_M(s-1)+i^{1, M}计算φ_Ms+i^{1, M}(s=2, 3, …, 2^M-1-1)。

2.4 Hilbert曲线初始条件表

在3维欧拉空间中，根据2.3节步骤 (1)、(2)、(3) 可以计算3维1阶Hilbert曲线的对应的地址顺序列表 (总共3×2^3-1=12条)。这12条Hilbert曲线是获取3维N阶Hilbert曲线经过网格点的地址顺序列表的初始条件，称为初始条件表。实际上，这12条Hilbert曲线就是构建3维N阶Hilbert曲线的基单元。

通过上面研究已经确定一个初始条件表，为了生成不同剖分等级N的Hilbert曲线，还需要研究两个相邻剖分等级之间的关系。

2.5 Hilbert曲线层级演进关系表

初始条件表确定地址顺序列表总体走向，而N阶和N+1阶之间的层级演进关系表可以确定子超立方体的分布情况，从而得到完整的Hilbert曲线地址顺序列表。N阶和N+1阶之间的层级演进关系是通过1阶和2阶之间的演进层级关系归纳总结而来。在M(M≥2) 维欧拉空间R^M中，层级演进关系表也描述了两个邻近超立方体的地址φ_l^{N, M}(h) 和φ_l^{N, M}(h+1) 的连接关系。

在子超立方体S_a^{N, M}中，采用地址序列{W_i^{n, M}}(n=1, 2, …, N-1 i=1, 2, …, N2^2N-1) 对应于第l条Hilbert曲线地址序列集合{φ_l^{n, M}}(l=1, 2, …, N2^N-1)。地址序列W_i^{n, M}包括起点和终点地址。N阶和N+1阶之间的演进层级关系可以使用W_i^{n, M}来描述，即

式中，l′=(l-1)2^M; l=1, 2, …, N2^N-1。如果W为空 (n=0)，则式 (8) 就等价于初始条件表。

{W_i^{n, M}}具有如下的性质

根据上述性质，可将演进层级关系公式简化为

式中，l″=(l-1)(2^M-1+1)，l=1, 2, …, N2^N-1。假定记演进层级关系表为T_induct^M[l]，则

式中，induct (W) 表示输入地址序列W，输出其对应第几条Hilbert曲线。

以3维欧拉空间R³为例，则T_induct^3[1]=(3, 2, 2, 10, 10, 5, 5, 9)，T_induct^3[3]=(1, 2, 2, 6, 6, 1, 1, 7)。根据T_induct³[l]可以计算得到其完整演进层级关系表达式，即3维Hilbert曲线N阶和N+1阶之间关系为

从式 (12) 可得到演进层次表T_induct³为

由此，可以根据初始条件表和演进层次表计算3维欧拉空间R³任意一条可能的Hilbert曲线经过的坐标地址顺序列表。当然，其也可以拓展到更高的维度空间。

为了与表 2对应，3维1阶Hilbert曲线的走向选定为第2条Hilbert曲线，即φ₂^{1, 3}为3维1阶Hilbert曲线对应的地址顺序列表，其中J对应X₁，I对应X₂，K对应X₃。3维N阶Hilbert曲线是在第2条Hilbert曲线基础上依次按照T_induct³得到的，其相应的地址顺序列表 (即Hilbert编码) 也可以计算得到。

3 基于Hilbert曲线层级演进关系的网格编码计算

全球立体网格可抽象为一个初始超立方体，记为S₀^{0, 3}；按照八叉树剖分N次得到的立体网格单元可以抽象为子超立方体，记为S_a^{N, 3}；φ^{N, 3}(h) 表示Hilbert编码为h的立体网格单元。若Hilbert编码以二进制形式表达，则记为h²；若以八进制形式表达，则记为h⁸。

八叉树立体网格Hilbert编码是一种可用{0, 1, 2, …, 7}8个码元描述的立体网格单元“定位”系统。Hilbert编码h²的每3位对应于立体网格单元S_a^{n, 3}的一次“八分”，从最高的3位到最低3位依次对应S_a^{0, 3}的第1级细分、第2级细分……直至第N级细分。

由于地球立体空间是一个连续的空间，因此可不断增加剖分等级N，减小立体网格单元的粒度，提高网格定位的精度。随着剖分等级N不断增加，八叉树立体网格Hilbert编码集合H将变成一个满足Hilbert曲线生成法则的无限空间集合。结合一般数学加法运算的性质，根据分析，Hilbert编码加法满足下面的性质：

(1) 封闭性：若h_A², h_B²∈H，h_new²=h_A²+h_B²，则h_new²∈H。

(2) 结合律：若h_A², h_B², h_c²∈H，则满足 (h_A²+h_B²)+h_c²=h_A²+(h_B²+h_c²)。

(3) 交换律：若h_A², h_B²∈H，则满足h_A²+h_B²=h_B²+h_A²。

(4) 消去律：若h_A², h_B², h_c²∈H且h_A²+h_c²=h_B²+h_c²，则h_A²=h_B²。

假定A和B是两个同一剖分等级的立体网格单元，分别为S_a^{N, 3}和S_a′^{N, 3}。根据Hilbert编码生成算法和层级关系，可以得出Hilbert编码满足下列运算。

3.1 集合运算

(1) 若a≠a′，则S_a^{N, 3}≠S_a′^{N, 3}或S_a^{N, 3}∩S_a′^{N, 3}=∅，即两个不同立体网格单元的Hilbert编码不同，地址也不同。

(2) ，S^{N+1, M}_abi是S_a^{N, 3}经过一次剖分得到8个子单元。

3.2 位运算

(1) 异或运算。h_A²⊕h_B²，表示A和B两个立体网格单元的Hilbert编码逐位进行异或运算。异或运算可以应用于判断Hilbert编码在哪一个剖分等级开始与其直接父单元不再是同一个单元。

(2) 右移运算。h_N+1²≫3=h_N²，即N+1级立体网格单元S^{N+1, 3}_abi(0≤i≤7) 的Hilbert编码右移3位得到其对应的N级单元S_a^{N, 3}的Hilbert编码。

(3) 左移运算。h_N²≪3=h_N+1²，即N级立体网格单元S_a^{N, 3}的Hilbert编码左移3位得到N+1级单元S^{N+1, 3}_ab0(第1个子单元) 的Hilbert编码。

3.3 数学加法运算

h_new²=h_A²+h_B²或h_new⁸=h_A⁸+h_B⁸，表示两个Hilbert码进行数学加法运算，其运算规则按照二进制或八进制加法进行。Hilbert编码 (h⁸) 的每1位表示立体网格单元在相应剖分层次上的8个子立方体中的哪一个。

具体运算过程 (以八进制为例)：从低位开始逐位进行加法运算，若相应结果大于等于8，则向高位进位，其含义是在高一等级的剖分网格中沿Hilbert曲线走向前进了相应的子立方体数量；若最终结果增加了一位，则表明剖分等级增加一级。若想得到Hilbert编码对应的地址码φ^{N, 3}(h)，只需依次从高位到低位获取数值，依据Hilbert编码层级关系计算得到。

例如：44+72=136，计算结果增加一位表明剖分等级增加一级，个位数“6”表示在最低剖分等级的8个子立方体的第7个中 (注意：编码是从0，1，…，7)，同理依次可以得出10位数“3”对应于第4个子立方体，百位数“1”对应于第1个子立方体。按照层次关系，可得知不同剖分层次上对应的Hilbert曲线依次为φ₂^{3, 3}、φ₁^{2, 3}、φ₂^{1, 3}，则可得出136(Hilbert码) 对应的地址码φ^{3, 3}(136)=123，二进制形式为001010011。

3.4 数学减法运算

h_new²=|h_A²-h_B²|或h_new⁸=|h_A⁸-h_B⁸|，表示两个Hilbert编码进行数学减法运算，其运算规则按照二进制或八进制减法进行。

具体运算过程 (八进制为例)：先比较两个Hilbert编码的大小，以较大的Hilbert编码减去较小的，从低位开始逐位进行减法运算，若不能进行减法运算，则向高位借位，其含义是在高一等级的剖分网格中沿Hilbert曲线走向后退一个子立方体；若最终结果减少了一位，则表明剖分等级减少一级。若想得到Hilbert编码对应的地址码φ^{N, 3}(h)，只需依次从高位到低位获取数值，依据Hilbert编码层级关系计算得到。

4 结束语

建立以网格编码为基础的代数运算法则是推动全球空间网格从数据组织管理框架发展成为新一代空间分析模型的理论基础。本文基于Hilbert曲线设计了八叉树立体网格单元编码，分析推导出Hilbert曲线层级演进关系表，以此基础给出了用于网格编码代数运算的若干操作算子。Hilbert曲线层级演进关系的分析方法可以扩展到其他结构的全球空间网格模型中，从而为构建全球空间网格分析理论与方法提供参考。