Hilbert曲线先将$n$维空间规则划分为无缝无叠的${\left({\mathrm{2}}^m\right)}^n$个超立方体网格单元（称为格元），再利用一维曲线不遗漏、不交叉地穿行所有格元，实现了$n$维空间$R^n$与一维空间$R^{\mathrm{1}}$的一一映射。三维Hilbert曲线实现了三维空间至一维空间的一一映射，其对边长为$T$的立方体空间作嵌套八叉剖分，第一次剖分得到8个子格元，经过$m$次剖分后得到${\mathrm{8}}^m$个子格元，每个子格元边长为$T/{\mathrm{2}}^m$，$m$层级格元对应于$m$阶曲线。层级$m$的取值范围为$\left[\mathrm{0},+\infty \right.)$，当$m=\mathrm{0}$时，格元为根格元。

将曲线在各层级格元中的填充顺序$\epsilon$相互衔接作为$m$层级格元的Hilbert码，表示为$h={\epsilon }^{\mathrm{1}}\oplus {\epsilon }^{\mathrm{2}}\oplus \dots \oplus {\epsilon }^m$，如图 1所示。因此，Hilbert码涵盖了两个重要层级剖分性质：（1）Hilbert码各个位上的数值为曲线在各层级格元中的填充顺序；（2）两个同属于一个父格元的子格元，其Hilbert码除最后一位外前缀完全相同。

图  1  一层级、二层级格元及其Hilbert码

Figure 1.  First and Second Level Grid Elements and Their Hilbert Codes

将层级为$m$、Hilbert码为$h$的格元记作$v_h^m$。根据上述曲线与格元性质，作以下计算、概念约定：（1）与格元$v_h^m$层级差为$k$的父格元记作$P^k\left(v\right)=v_{h\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{8}}^{m-k}$，其中$h\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{8}=⌊h/{\mathrm{8}}^k⌋$，表示除去Hilbert码$h$的最后$k$位数。本文中的父格元指与格元$v_h^m$层级差为$\mathrm{1}$的父格元，简记为$P\left(v\right)$。（2）格元$v_h^m$进行一次剖分得到的8个子格元中，填充顺序为$i$的子格元$C\left(v,i\right)=v_{h\oplus i}^{m+\mathrm{1}}$，$C\left(v,i\right)$称作格元$v$的第$i$个子格元。（3）格元$v_h^m$在$P\left(v\right)$剖分得到的8个子格元中的填充顺序$I\left(v\right)=h\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{8}$，其中$\mathrm{m}\mathrm{o}\mathrm{d}$为取余运算，$I\left(v\right)$即为Hilbert码$h$的最后1位数。（4）为区分不同的邻近关系，本文将同一层级格元间的邻近关系分为兄弟邻近以及堂兄弟邻近两种关系。由同一格元经过一次剖分得到的称为兄弟邻近格元，其余称为堂兄弟邻近格元。例如格元$v_{\mathrm{30}}^{\mathrm{2}}$与$v_{\mathrm{31}}^{\mathrm{2}}$、$v_{\mathrm{32}}^{\mathrm{2}}$、$v_{\mathrm{36}}^{\mathrm{2}}$为兄弟邻近关系，与$v_{\mathrm{47}}^{\mathrm{2}}$为堂兄弟邻近关系。兄弟邻近格元Hilbert码除最后一位外前缀完全相同。

为区分不同的空间方向D，将空间方向共分为上（top，T）、下（bottom，B）、左（left，L）、右（right，R）、内（in，I）、外（out，O）6个方向。格元进行一次剖分得到8个子格元，将其按照方向位置分别称为LOB、LOT、ROT、ROB、RIB、RIT、LIT、LIB，如图 2所示。线性八叉树邻近格元分为面相邻的面邻近、边相邻的边邻近、角相邻的角邻近3种，例如LOT为ROT的面邻近，LOB为ROT的边邻近，LOT为RIB的角邻近。在同层级下，线性八叉树任一格元的邻近格元集合是由6个面邻近、12个边邻近、8个角邻近等共26个邻近格元组成。本文主要研究的邻域计算问题就是快速计算任一格元的26个邻近格元Hilbert码。

图  2  一次剖分得到的8个子格元

Figure 2.  8 Sub-Grid Elements Obtained by One Split

本文的研究目的是仅使用Hilbert曲线自身性质实现邻近格元Hilbert码计算。因此，了解掌握Hilbert曲线的填充规律是进行邻近Hilbert码计算的前提。首先对Hilbert曲线的填充规律进行分析，提出状态矩阵以及层级演进退化函数概念，然后将其应用于邻近格元Hilbert码的计算中。

Hilbert曲线的自相似与自复制特征决定了Hilbert曲线中只有24种形态的基元曲线，不管多复杂的Hilbert曲线都由这24种基元曲线递归派生而来^[22]，分别为${\phi }_{\mathrm{1}},{\phi }_{\mathrm{2}}\cdots {\phi }_{\mathrm{24}}$，如图 3所示，每种形态基元曲线的区分标准是Hilbert曲线在由一个格元剖分得到的LOB、LOT、ROT、ROB、RIB、RIT、LIT、LIB中的填充状态。

图  3  Hilbert基元曲线

Figure 3.  Hilbert Primitive Curve

为描述这24种形态基元曲线，引入状态矩阵S记录24种基元曲线，按照LOB、LOT、ROT、ROB、RIB、RIT、LIT、LIB的顺序，记录第$k$种基元曲线在填充至当前子格元时共经过了多少个子格元，记作状态向量${\mathit{s}}_k\left(k\in \left(\mathrm{1,\; 2}\cdots \mathrm{24}\right)\right)$，然后将24个行向量合成矩阵S，矩阵S的第$u$行对应状态向量${\mathit{s}}_{u+\mathrm{1}}$，对应基元曲线为${\phi }_{u+\mathrm{1}}$。

针对一个格元，Hilbert曲线在其子格元间的填充状态为${\phi }_k$，Hilbert曲线在各个子格元再次剖分后得到的格元间的填充状态由${\phi }_k$决定，这种映射关系不随层级变化而发生改变^[23]。

如图 4所示，二阶曲线的填充状态为${\phi }_{\mathrm{2}}$，对应状态向量为s₂；剖分后Hilbert曲线在各个子格元间的填充状态分别为${\phi }_{\mathrm{3}}\mathrm{、}{\phi }_{\mathrm{1}}\mathrm{、}{\phi }_{\mathrm{2}}\mathrm{、}$${\phi }_{\mathrm{10}}\mathrm{、}{\phi }_{\mathrm{3}}\mathrm{、}{\phi }_{\mathrm{2}}\mathrm{、}{\phi }_{\mathrm{5}}\mathrm{、}{\phi }_{\mathrm{10}}$，状态向量分别对应${\mathit{s}}_{\mathrm{3}}\mathrm{、}{\mathit{s}}_{\mathrm{1}}\mathrm{、}{\mathit{s}}_{\mathrm{2}}\mathrm{、}{\mathit{s}}_{\mathrm{10}}\mathrm{、}{\mathit{s}}_{\mathrm{3}}\mathrm{、}{\mathit{s}}_{\mathrm{2}}\mathrm{、}{\mathit{s}}_{\mathrm{5}}\mathrm{、}{\mathit{s}}_{\mathrm{10}}$，则存在状态向量层级映射关系$T({\mathit{s}}_{\mathrm{2}})\to$${\mathit{s}}_{\mathrm{3}},{\mathit{s}}_{\mathrm{1}},$ ${\mathit{s}}_{\mathrm{2}},{\mathit{s}}_{\mathrm{10}},{\mathit{s}}_{\mathrm{3}},{\mathit{s}}_{\mathrm{2}},{\mathit{s}}_{\mathrm{5}},{\mathit{s}}_{\mathrm{10}}$。这种状态向量映射关系并不会随层级变化而发生改变，即无论格元层级的高低，若其状态向量为${\mathit{s}}_{\mathrm{2}}$，其子格元的状态向量均为${\mathit{s}}_{\mathrm{3}}\mathrm{、}{\mathit{s}}_{\mathrm{1}}\mathrm{、}{\mathit{s}}_{\mathrm{2}}\mathrm{、}{\mathit{s}}_{\mathrm{10}}\mathrm{、}{\mathit{s}}_{\mathrm{3}}\mathrm{、}{\mathit{s}}_{\mathrm{2}}\mathrm{、}{\mathit{s}}_{\mathrm{5}}\mathrm{、}{\mathit{s}}_{\mathrm{10}}$。

图  4  状态向量演进

Figure 4.  State Vector Evolution

引入状态演进矩阵E记录24种状态向量与下一层级格元状态向量的映射关系。E第u行$\left(u\in \left(\mathrm{0,\; 1},\mathrm{2}\cdots \mathrm{23}\right)\right)$对应映射关系$T\left({\mathit{s}}_{u+\mathrm{1}}\right)$，即：

根据E可知，若已知格元$v_h^m$的状态向量${\mathit{s}}_k$，可求得其填充顺序为$i$的子格元的状态向量$\mathit{s}\left(C\left(v,i\right)\right)$为：

式中，$T^c\left({\mathit{s}}_k,i\right)$为状态演进函数。使用状态演进函数对格元$v_h^m$进行一次演进，可得到格元的第$i$个子格元的Hilbert码以及状态向量。借助函数$T^c\left({\mathit{s}}_k,i\right)$，以格元$v_h^m$的Hilbert码$h$为输入，计算格元$v_h^m$状态向量的算法流程如图 5所示，其中$h_{\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{g}\mathrm{t}\mathrm{h}}$为Hilbert码的长度，$h_{\mathrm{b}\mathrm{i}\mathrm{t}}\left(i\right)$为Hilbert码$h$的第$i$位。

图  5  状态向量计算流程

Figure 5.  Process of State Vector Calculation

对函数$T^c\left({\mathit{s}}_k,i\right)$求逆，可得其逆函数为：

式中，$\mathit{s}\left(C\left(v,i\right)\right)$为格元$v_h^m$的第$i$个子格元的状态向量；$T^p\left(\mathit{s}\left(C\left(v,i\right)\right),i\right)$为状态退化函数。使用状态退化函数对子格元$C\left(v,i\right)$进行一次退化，可得到格元$v_h^m$的Hilbert码以及状态向量。

§1.1中已经按照空间方向对一次剖分得到的8个子格元的位置进行区分，由此推导得到邻近方向D上各个子格元的同层级面邻近子格元位置，结果如表 1所示

因此，根据Hilbert曲线所涵盖的层级剖分性质，可在已知中心格元状态向量的前提下，计算面邻近格元Hilbert码，步骤如下：（1）以中心格元为起点，层级退化搜索共同父格元。（2）以共同父格元为起点，层级演进计算得到面邻近格元的Hilbert码。层级退化演进过程中，利用状态向量判断当前格元所处的位置，而层级退化演进函数则用于计算不同层级格元的状态向量。每一次退化都将中心格元在D方向上的邻近格元位置存入数组Direction[ ]中。演进时，每一次选择的演进格元位置只需从Direction[ ]的末尾开始往前取值即可。此外，还需在各次层级退化演进中检查八叉树中是否存在当前格元。面邻近格元Hilbert码计算算法具体流程如图 6所示。

图  6  面邻近格元Hilbert码计算流程

Figure 6.  Calculation Process of the Hilbert Code of the Neighboring Grid Element

以计算中心格元$v_{\mathrm{042}}^{\mathrm{3}}$在R方向上面邻近格元Hilbert码为例。详细阐述本文面邻近格元Hilbert码计算过程，步骤如下：（1）以格元$v_{\mathrm{042}}^{\mathrm{3}}$为起点退化一层级，计算得到其父格元$v_{\mathrm{04}}^{\mathrm{2}}$的状态向量为${\mathit{s}}_{\mathrm{2}}$，并推得$v_{\mathrm{042}}^{\mathrm{3}}$属于ROT。查询表 1可知，ROT在R方向上的面邻近格元为LOT*，LOT*与ROT不属于同一个父格元，LOT*与ROT之间属于堂兄弟面邻近关系，无法使用格元$v_{\mathrm{04}}^{\mathrm{2}}$演进一层级的方法计算LOT*的Hilbert码。（2）继续以格元$v_{\mathrm{04}}^{\mathrm{2}}$为起点退化一层级，计算得到其父格元$v_{\mathrm{0}}^{\mathrm{1}}$的状态向量为${\mathit{s}}_{\mathrm{3}}$。在状态向量${\mathit{s}}_{\mathrm{3}}$中，填充顺序为$\mathrm{4}\left(\mathrm{0}\underline{\mathrm{4}}\mathrm{2}\right)$的子格元为ROT，$v_{\mathrm{04}}^{\mathrm{2}}$在R方向上的面邻近格元仍然与$v_{\mathrm{04}}^{\mathrm{2}}$不同属于一个父单元。（3）继续以格元$v_{\mathrm{0}}^{\mathrm{1}}$为起点退化1层级，计算得到根格元$v^{\mathrm{0}}$的状态向量为${\mathit{s}}_{\mathrm{2}}$。在状态向量${\mathit{s}}_{\mathrm{2}}$中，填充顺序为$\mathrm{0}\left(\underline{\mathrm{0}}\mathrm{42}\right)$的子格元为LOB，LOB在R方向上的子格元是ROB，LOB与ROB同属于一个父格元，LOB与ROB之间属于兄弟面邻近关系，可使用演进1层级的方法计算LOB的Hilbert码。（4）在状态向量${\mathit{s}}_{\mathrm{2}}$中，ROB的填充顺序为3。以根格元$v^{\mathrm{0}}$为起点，向下演进1层级，得到根格元$v^{\mathrm{0}}$填充顺序为3的子格元$C\left(v^{\mathrm{0}},\mathrm{3}\right)=v_{\mathrm{3}}^{\mathrm{1}}$，$v_{\mathrm{3}}^{\mathrm{1}}$的状态向量$\mathit{s}\left(v_{\mathrm{3}}^{\mathrm{1}}\right)=\mathit{s}\left(C(v^{\mathrm{0}},\mathrm{3})\right)=s_{T^c\left({\mathit{s}}_{\mathrm{2}},\mathrm{3}\right)}=s_{\mathit{E}[\mathrm{1}][\mathrm{3}]}={\mathit{s}}_{\mathrm{10}}$。（5）在状态向量${\mathit{s}}_{\mathrm{10}}$中，LOT的填充顺序为1。以格元$v_{\mathrm{3}}^{\mathrm{1}}$为起点，向下演进1层级，计算格元$v_{\mathrm{3}}^{\mathrm{1}}$填充顺序为1的子格元$C\left(v_{\mathrm{3}}^{\mathrm{1}},\mathrm{1}\right)=v_{\mathrm{31}}^{\mathrm{2}}$，$v_{\mathrm{31}}^{\mathrm{2}}$的状态向量$\mathit{s}\left(v_{\mathrm{31}}^{\mathrm{2}}\right)=\mathit{s}\left(C(v_{\mathrm{3}}^{\mathrm{1}},\mathrm{1})\right)={\mathit{s}}_{\mathrm{13}}$。（6）在状态向量${\mathit{s}}_{\mathrm{13}}$中，LOT的填充顺序为1。以格元$v_{\mathrm{31}}^{\mathrm{2}}$为起点，向下演进一层级，计算格元$v_{\mathrm{31}}^{\mathrm{2}}$填充顺序为1的子格元$C\left(v_{\mathrm{31}}^{\mathrm{2}},\mathrm{1}\right)=v_{\mathrm{311}}^{\mathrm{3}}$，$v_{\mathrm{311}}^{\mathrm{3}}$的状态向量$s\left(v_{\mathrm{311}}^{\mathrm{3}}\right)=\mathit{s}\left(C(v_{\mathrm{31}}^{\mathrm{2}},\mathrm{1})\right)={\mathit{s}}_{\mathrm{10}}$。检测八叉树中不存在格元$v_{\mathrm{311}}^{\mathrm{3}}$，则停止演进，返回Hilbert码$\mathrm{31}$，最终计算得到格元$v_{\mathrm{042}}^{\mathrm{3}}$在R方向上面邻近格元为$v_{\mathrm{31}}^{\mathrm{2}}$。整体计算过程如图 7所示，黄色格元为中心格元$v_{\mathrm{042}}^{\mathrm{3}}$。

图  7  计算R方向上面邻近格元Hilbert码

Figure 7.  Hilbert Code Calculation of Adjacent Grid Elements in the R Direction

分析上述面邻近格元Hilbert码计算算法流程可知，面邻近格元Hilbert码计算算法的演进与退化过程都最多为$k$次循环（$k$为中心格元$v_h^m$与共同父格元之间的层级差），因此计算同层级邻近格元Hilbert码时，面邻近格元Hilbert码计算算法的时间复杂度为$O\left(\mathrm{2}k\right)$。

相比现有的转换算法^[13]，本文算法的优势如下：（1）定义了状态向量与层级演进退化函数，算法中不需要将Hilbert码转换为其他形式的坐标或空间排列码，避免了复杂Hilbert码层级拆解转换运算，不存在多种分支判断计算或嵌套循环计算，算法计算复杂度较小。（2）算法只具有单层循环，共需要进行两次步长为$k$的循环计算，算法的时间复杂度为$O\left(\mathrm{2}k\right)$。相比于转换算法的时间复杂度$O\left(\mathrm{6}m\right)$，层级差$k$必定不大于中心格元层级$m$，故本文算法的循环计算次数必定不大于转换算法，且本文算法的计算复杂度更小，整体算法效率更高。

角邻近格元Hilbert码计算过程可以由3次面邻近格元Hilbert码计算组成，如图 8所示，首先以黄色格元为中心格元，计算B方向上的面邻近无色格元Hilbert码，然后再以无色格元为中心格元计算L方向上的面邻近灰色格元Hilbert码，最后再以灰色格元为中心格元计算I方向上的面邻近绿色格元Hilbert码，得到黄色格元角邻近格元中的绿色格元Hilbert码。

图  8  角邻近格元Hilbert码计算

Figure 8.  Hilbert Code Calculation of Corner Adjacent Grid Element

为了验证算法正确性、性能与效率，设计了邻近格元计算正确性、状态向量计算效率、Hilbert码计算效率3组实验进行测试。实验硬件环境为CPU Intel Core i7-7700K（双核4.2 GHz），内存为64 GB，硬盘速度为7 200 r/min；软件环境为Visual Studio 2015，Release版本，X64，C++。

Hilbert曲线具有空间填充曲线的分形特征，具有自相似与自复制特性，高层级曲线由低层级曲线平移、求反与交换推导得出。对比不同算法的计算结果，若低层级结果一致，由分形特征可知，其高层级结果也一致^[24]。因此，使用本文算法与文献[13]算法分别计算不同层级中心格元邻近格元的Hilbert码并进行对比验证。

将层级为$m$的所有网格单元集合记为$V^m$，将文献[13]算法求得的Hilbert码记作$I(v)$，将本文算法生成的Hilbert码记作$H(v)$。对于同一邻近格元，定义两种算法计算Hilbert码差值$\mathrm{\Delta }\mathrm{C}\mathrm{o}\mathrm{d}\mathrm{e}(v)$为：

定义$m$层级单元的Hilbert码差值之和$G\left(V^m\right)$为：

计算了1~10层的Hilbert码差值之和，结果如表 2所示。由表 2可知，各层级$G\left(V^m\right)$均为0，证明本文算法与文献[13]算法生成的Hilbert码结果一致。10层级之内本文算法结果正确，由分形特征推知更高层级结果仍正确。

相比转换算法，本文算法执行之前需要已知中心格元的状态向量。为探究中心格元状态向量计算步骤的执行效率是否会对后续邻近格元Hilbert码计算产生影响，设计实验进行测试。在1~20层级的所有格元中，随机挑选$\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$个格元（存在重复格元）进行状态向量计算，记录其计算总耗时，并计算各个层级内的计算效率，即1 ms时间内可进行计算的次数。图 9为1~20层级状态向量计算耗时与计算效率，展现了计算总时间与计算效率随层级变化的趋势。

图  9  状态向量计算耗时与计算效率

Figure 9.  Time-Consuming and Computational Efficiency of State Vector Calculation

由图 9可知，计算平均时间对应的曲线随着层级的增加而逐步上升，表明进行一次状态向量计算消耗的时间与层级呈正相关，其中第20层级进行一次计算的时间为$\mathrm{238}\mathrm{\mu }\mathrm{s}$。计算效率对应曲线随层级的增加而逐步下降，表明状态向量计算效率与层级呈负相关，其中第20层级1 ms内仍可对4 201个格元进行计算。

为分析本文算法计算邻近格元Hilbert码的效率性能，将本文算法与现有的Morton码转换算法^[13]进行对比实验。实验1：对比相同层级下，两种算法计算不同次数所需的总耗时；实验2：对比相同计算次数下，两种算法在不同层级下计算所需的总耗时。

在层级15情况下，随机选取当前层级的$n$个格元，计算其状态序列，并采用两种算法计算26个邻近格元Hilbert码。实验过程中，$n$的取值分别为$\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}\cdots \mathrm{8}\times {\mathrm{10}}^{\mathrm{6}}$，统计各数量格元邻近格元Hilbert码的总耗时，结果如表 3所示，查询时间随格元个数变化如图 10所示。由图 10可知，两种算法对应曲线随着格元个数的增加均呈现近似线性上升，表明两种算法的计算所需时间均与格元个数呈正相关；对比图 10中相同格元数下两种算法所需的计算时间，本文算法的计算时间均要少于转换算法，计算速度为转换算法的2.1~2.4倍。

图  10  查询时间随格元个数变化

Figure 10.  Query Time Varies with the Number of Grid Elements

在层级1~20下，随机选取当前层级中的$\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$个格元作为样本集，计算样本集内所有格元状态序列，并分别采用两种算法计算其26个邻近格元的Hilbert码，统计两种算法计算Hilbert码总耗时后，对一次邻域计算的平均时间进行计算，结果如图 11所示。将Morton码转换算法的计算时间除以本文算法的计算时间，作为本文算法效率提升的倍数，各层级结果如表 4所示。

图  11  计算时间对比

Figure 11.  Comparison of Calculation Time

由图 11可知，转换算法对应曲线随层级增加呈现较为均匀的增长趋势，表明转换算法的邻近格元Hilbert码计算效率随层级的增加而线性降低；本文算法对应曲线较为平缓，随着层级的增加存在较为微小的增长，但基本处于100 $\mathrm{\mu }\mathrm{s}$以下，表明本文算法受层级影响较小，更加满足线性八叉树邻域计算需求。对比图 11中相同层级下两种算法所需的计算时间，本文算法所需时间均小于转换算法，且本文算法较转换算法的效率提升倍数随层级增加有扩大趋势，由表 4可知，在第20层级，本文算法效率相比于转换算法提升了2.6倍。

由§3.2中实验结果可知，计算一个格元的状态向量所需时间的数量级为${\mathrm{10}}^{-\mathrm{4}}\mathrm{\mu }\mathrm{s}$，进而可推得26个格元的状态向量所需时间的数量级为${\mathrm{10}}^{-\mathrm{3}}\mathrm{\mu }\mathrm{s}$，而进行一次邻近Hilbert计算所需平均时间数量级为${\mathrm{10}}^{\mathrm{2}}\mathrm{\mu }\mathrm{s}$，对比两者可知，计算格元状态向量所需时间可在计算邻近Hilbert码时忽略不计。

综上所述，在计算邻近格元Hilbert码的效率上，本文算法的效率要优于转换算法。

为提高基于Hilbert曲线编码的线性八叉树邻近格元计算效率，本文基于Hilbert曲线的填充规律设计了一种邻近格元Hilbert码的快速计算算法，并详细阐述了该算法的流程。两组实验对比结果表明，本文算法在邻近格元Hilbert码计算的效率上优于现有的Morton码转换算法，可作为高效数据检索、三维拓扑关系确定与连通性判断等空间操作的基础。