引用本文
贾连印, 孔明, 王维晨, 李孟娟, 游进国, 丁家满. 数据偏斜分布下的二维Hilbert编解码算法[J]. 清华大学学报(自然科学版), 2022, 62(9): 1426-1434.  
JIA Lianyin, KONG Ming, WANG Weichen, LI Mengjuan, YOU Jinguo, DING Jiaman. 2-D Hilbert encoding and decoding algorithms on skewed data[J]. Journal of Tsinghua University (Science and Technology), 2022, 62(9): 1426-1434.  
数据偏斜分布下的二维Hilbert编解码算法
贾连印1,3, 孔明1, 王维晨1, 李孟娟2, 游进国1,3, 丁家满1,3    
1. 昆明理工大学 信息工程与自动化学院, 昆明 650500;
2. 云南师范大学 图书馆, 昆明 650500;
3. 昆明理工大学 云南省人工智能重点实验室, 昆明 650500
收稿日期：2021-06-28
基金项目：国家自然科学基金地区基金项目(62162038, 62062046, 61562054)
作者简介：贾连印(1978－), 男, 副教授
通讯作者：丁家满, 副教授, E-mail：jiamanding@kust.edu.cn
摘要：高效的Hilbert曲线的编解码算法作为Hilbert曲线应用的基础，具有重要的研究意义。现有算法多未考虑数据偏斜分布的影响，因此在数据偏斜分布时效率较低。该文发现：对于特定的前m阶坐标，其对应的前m阶编码值与其第1阶编码值呈现特定的倍数关系；对于特定的前m阶编码值，其对应的前m阶坐标与其第1阶坐标呈现特定的倍数关系。基于这一发现，在融合高效位操作、快速置位检测等技术的基础上，提出了跳过前m阶的编码(skipping the first m orders Hilbert encoding, SFO-HE)算法和跳过前m阶的解码(skipping the first m orders Hilbert decoding, SFO-HD) 算法。这2个算法无需对前m阶逐阶编解码，可有效提高数据向Hilbert空间4个顶点偏斜时的编解码效率。扩展实验表明：该文算法对数据偏斜分布具有更好的适应性，在特定偏斜分布时效率大幅优于现有算法。
关键词：Hilbert曲线    状态视图    偏斜分布    编解码算法    
2-D Hilbert encoding and decoding algorithms on skewed data
JIA Lianyin1,3, KONG Ming1, WANG Weichen1, LI Mengjuan2, YOU Jinguo1,3, DING Jiaman1,3    
1. Faculty of Information Engineering and Automation, Kunming University of Science and Technology, Kunming 650500, China;
2. Library, Yunnan Normal University, Kunming 650500, China;
3. Yunnan Key Laboratory of Artificial Intelligence, Kunming University of Science and Technology, Kunming 650500, China
Abstract: Hilbert encoding and decoding are fundamental steps in many Hilbert curve based applications. However, existing algorithms are not very efficient when the data distribution is skewed. This paper shows that for a coordinate with the specific first m orders, the code of the first m orders is a multiple of its corresponding first order code. For a code with the specific first m orders, the coordinate of the first m orders is a multiple of its corresponding first order coordinate. These findings were used to develop an algorithm that skips the first m orders of the Hilbert encoding (SFO-HE) and another algorithm that skips the first m orders of the Hilbert decoding (SFO-HD). These algorithms exploit efficient bit operations and fast bit set detections to improve the encoding and decoding efficiencies for data skewed to the 4 corners of the Hilbert space. Extensive tests show that these two algorithms have good skewness adaptability and outperform existing algorithms on specific skewed data.
Key words: Hilbert curve    state view    skewed distribution    encoding and decoding algorithms    

空间填充曲线可将多维数据映射为一维数据，同时保持其在多维空间的局部性，使得复杂的多维空间问题可在一维空间下高效解决，因此在空间数据处理[1-2]、图像处理[3-4]、数据挖掘[5]等诸多领域得到了广泛的应用。

Z曲线[6]和Hilbert曲线[7]是最常见的2种空间填充曲线。Z曲线映射规则和编码简单，但跳变性大。Hilbert曲线的映射规则较复杂，但具有更好的空间连续性。高效的Hilbert曲线的编解码算法是Hilbert应用的基础，具有重要研究价值。本文着重于二维Hilbert曲线的编解码，因其应用最广泛。

二维Hilbert编解码算法研究由来已久，许多高效的算法[8-10]被提出。经典的二维Hilbert编码算法可划分为递归和迭代两类。文[8]设计了一种以字节为导向的Hilbert曲线生成算法，但由于是基于递归的方法，需要较大的映射开销，效率并不理想。文[11]指出多数递归算法效率较低。因此后续研究多基于迭代算法来展开。

迭代算法中，基于状态视图的编码算法由于便捷高效、易实现等优点得到了广泛深入的研究。文[12]提出了复杂度为O(n)的基于状态矩阵的二维Hilbert快速编码算法，构建了高效的状态视图，并通过迭代查询状态视图进行编码。文[13]通过象限映射和逆映射，设计了一个统一状态视图，可同时支持编码和解码操作。

上述算法的主要不足在于：其对所有输入数据同等对待，未考虑偏斜分布的影响。为此，文[14]通过设计高效的状态视图并结合快速置位检测算法，提出了免计前0的二维Hilbert编码算法，从而提高了数据向原点偏斜时的编码效率。

此外，也有一些不依赖状态视图的算法实现。文[9]通过计算代替查找表，可降低状态视图的开销，但引入了较多计算开销。文[10]设计了基于位运算的高效非递归编码算法，能很好地支持二维及高维数据的编码。这2个算法被认为是目前Hilbert编码事实上的标准，得到了广泛的应用。

尽管文[14]的算法可在一定程度上解决数据偏斜分布的问题，但仅在数据向原点偏斜时高效。对于其他偏斜的场景，该算法不仅无法跳过阶，反而需要付出额外的置位检测代价，从而降低算法效率。

本文发现：对于特定的前m阶坐标，其对应的前m阶编码值与其第1阶编码值呈现特定的倍数关系；对于特定的前m阶编码值，其对应的前m阶坐标与其第1阶坐标呈现特定的倍数关系。基于这一发现，在融合高效位操作、快速置位检测等技术的基础上，提出了跳过输入前m阶的编码(skipping the first m orders Hilbert encoding, SFO-HE) 算法和跳过前m阶的解码(skipping the first m orders Hilbert decoding, SFO-HD) 算法。其中m对编码而言，表示输入坐标前部为连续0或1的阶数；对解码而言，表示输入编码值前部为连续00、01、10或11的阶数。对这m阶无需逐阶编解码，可有效提高数据特定偏斜分布下的编解码效率。

1 必要前提 1.1 Hilbert曲线的概念

Hilbert曲线对空间进行网格划分，1阶二维Hilbert曲线将整个空间分割成2×2共4个子区域，将左下、左上、右上、右下的4个子区域依次连接即可构成1阶开口朝下的二维Hilbert曲线。通过在每个网格嵌入一个1阶曲线，并对左下网格的曲线左旋90°并垂直翻转，对右下网格的曲线右旋90°并垂直翻转，即可得到2阶二维Hilbert曲线。采用类似的变换规则可进一步得到更高阶的二维Hilbert曲线。

1.2 符号表

本文涉及的符号如表 1所示。为便于表述，文中X=(X1X2…Xn)2表示X1X2…Xn为X的二进制形式，以此类推。

表 1 命名规则表
符号 描述
P(X, Y) 坐标分量X=(X1X2…Xn)2
Y=(Y1Y2…Yn)2
Z=(Z1Z2…Z2n－1Z2n)2 Hilbert编码值
n 阶数
m 跳过的阶数
Xi X的第i阶
Yi Y的第i阶
Xi, Yi 第i阶坐标值
Z2i－1Z2i Z的第i阶的编码值
Si 第i阶的状态
$\overline{X_{i}}=X_{1} X_{2} \cdots X_{i}$ X的前i阶
$\overline{Y_{i}}=Y_{1} Y_{2} \cdots Y_{i}$ Y的前i阶
$\overline{Z_{2 i-1} Z_{2 i}}=Z_{1} Z_{2} \cdots Z_{2 i-1} Z_{2 i}$ Z的前i阶的编码值
表选项

2 二维Hilbert曲线编码 2.1 状态视图的构建

状态视图是本文算法的重要基础，反映了物理坐标、二维Hilbert编码值和状态之间的联系。

本文状态视图基于1阶二维Hilbert曲线构建，不同的1阶曲线开口方向对应不同的1阶状态。根据开口方向不同，1阶二维Hilbert曲线共包含图 1的4种基本状态：状态0(开口向下)、状态1(开口向左)、状态2(开口向上)、状态3(开口向右)。图中括号内的数值为1阶物理坐标，物理坐标上的数值为对应的1阶编码值。

图 1 1阶二维Hilbert曲线对应的4种状态图
图选项

根据图 1，易得出状态Si下1阶物理坐标(Xi, Yi)与对应的编码值Z2i－1Z2i及下一阶状态Si+1之间的关系，并构建编码状态视图如表 2所示。

表 2 编码状态视图表
Si XiYi Z2i-1Z2i Si+1
0 00 00 1
0 01 01 0
0 10 11 3
0 11 10 0
1 00 00 0
1 01 11 2
1 10 01 1
1 11 10 1
2 00 10 2
2 01 11 1
2 10 01 2
2 11 00 3
3 00 10 3
3 01 01 3
3 10 11 0
3 11 00 2
表选项

为降低空间开销并提高编码效率，将表 2设计为2个3维的数组：PZMap和PSMap。PZMap存储坐标值到编码值的映射，第i阶编码值Z2i－1Z2i可通过PZMap[Si][Xi][Yi]直接获取。PSMap则存储坐标值到下一阶状态的映射，第(i+1)阶状态Si+1可通过PSMap[Si][Xi][Yi]获取。为便于分析和实现，本文规定，编码前的初始状态S1=0。

2.2 编码 2.2.1 动机

一种比较直观的编码方式(如文[12]算法)是逐阶迭代检索创建的编码状态视图，以获得该阶编码值和下一阶的状态，进而获得最终的编码值。文[12]算法的复杂度为O(n)，其主要不足在于需要迭代访问所有的阶，算法效率较低。

为此，本文研究跳过输入数据前部全0或全1阶的方法，在提高编码效率的同时，使算法能适应数据偏斜到任意顶点的情形，提高算法的健壮性。

2.2.2 SFO-HE算法

SFO-HE算法旨在跳过输入数据前部全0或全1的m阶，基于以下基本定理。

定理1  若P(X, Y)的前m阶坐标($\overline{X_{m}}, \overline{Y_{m}}$)分别为(00…0, 00…0)2、(00…0, 11…1)2、(11…1, 00…0)2、(11…1, 11…1)2时，设定A=(4m－1)/3，则其前m阶编码值$\overline{Z_{2 m-1} Z_{2 m}}=\text { PZMap } \left[S_{1}\right]\left[X_{1}\right]\left[Y_{1}\right] \times A$，对应的第(m+1)阶的状态Sm+1分别为0(m为偶数时)或1(m为奇数时)、0、0(m为偶数时)或3(m为奇数时)、0。

证明：分4种情况分别讨论

1) ($\overline{X_{m}}, \overline{Y_{m}}$)为(00…0, 00…0)2的情形。

m=1时，由X1=0, Y1=0, S1=0，根据PZMap可知$\overline{Z_{1} Z_{2}}$=Z1Z2=(00)2=0, S2=1，定理成立。

m=2时，由X2=0, Y2=0, S2=1，根据PZMap可知Z3Z4=(00)2, S3=0，故$\overline{Z_{3} Z_{4}}=\overline{Z_{1} Z_{2}} Z_{3} Z_{4}=$(0000)2 =0，定理成立。

m≥3时，第3阶的输入坐标和状态与第1阶相同，故其输出编码和状态与第1阶相同。同理第4阶的输入和输出与第2阶相同，进入循环状态。故Z2m－1Z2m=(00)2，Sm+1=0(m为偶数)或1(m为奇数)，，故$\overline{Z_{2 m-1} Z_{2 m}}=$ PZMap $\left[S_{1}\right]\left[X_{1}\right]\left[Y_{1}\right] \times A$，定理成立。

2) ($\overline{X_{m}}, \overline{Y_{m}}$)为(00…0, 11…1)2的情形。

m=1时，由X1=0, Y1=1, S1=0, 根据PZMap可知$\overline{Z_{1} Z_{2}}$=Z1Z2=(01)2=1, S2=0，定理成立。

当m≥2时，可见每阶的输入坐标和状态均与第1阶相同，则输出编码和下一阶状态也均与第1阶相同，进入循环状态。故由Xm=0, Ym=1, Sm=0，易见Z2m－1Z2m=(01)2, Sm+1=0, $\overline{Z_{2 m-1} Z_{2 m}}=\overline{Z_{2 m-3} Z_{2 m-2}} Z_{2 m-1} Z_{2 m}=$ ，故$\overline{Z_{2 m-1} Z_{2 m}}={\rm{PZMap}}\left[S_{1}\right]\left[X_{1}\right]\left[Y_{1}\right] \times A$，定理成立。

3) ($\overline{X_{m}}, \overline{Y_{m}}$)为(11…1, 00…0)2的情形。

m=1时，由X1=1, Y1=0, S1=0，根据PZMap可知$\overline{Z_{1} Z_{2}}$=Z1Z2=(11)2=3, S2=3，定理成立。

m=2时，由X2=1, Y2=0, S2=3，根据PZMap可知Z3Z4=(11)2, S3=0，故$\overline{Z_{3} Z_{4}}=\overline{Z_{1} Z_{2}} Z_{3} Z_{4}=(1111)_{2}=15$，定理成立。

m≥3时，第3阶的输入坐标和状态与第1阶相同，故输出编码和状态与第1阶相同，同理第4阶的输入和输出与第2阶相同，进入循环状态。故Z2m－1Z2m=(11)2, Sm+1=0 (m为偶数)或3(m为奇数)，，故$\overline{Z_{2 m-1} Z_{2 m}}=\operatorname{PZMap}\left[S_{1}\right]\left[X_{1}\right]\left[Y_{1}\right] \times A$，定理成立。

4) ($\overline{X_{m}}, \overline{Y_{m}}$)为(11…1, 11…1)2的情形。

m=1时，由X1=1, Y1=1, S1=0，根据PZMap可知$\overline{Z_{1} Z_{2}}$=Z1Z2=(10)2=2, S2=0，定理成立。

m≥2时，可见每阶的输入坐标和状态均与第1阶相同，则输出编码和下一阶状态也均与第1阶相同，进入循环状态。故可得Z2m－1Z2m=(10)2, Sm+1=0, ，故$\overline{Z_{2 m-1} Z_{2 m}}=\mathrm{PZMap}\left[S_{1}\right]\left[X_{1}\right]\left[Y_{1}\right] \times A$，定理成立。

故定理得证。

由定理1可知，若输入数据前m阶为连续的0或1，则这m阶的编码及第(m+1)阶的状态可直接计算获得，无需逐阶编码，因此具有较高的编码效率。

快速检测输入数据前部为连续的0或1的阶数，即获取最大的m，对提高编码效率至关重要。高效的置位检测算法FFB (find first bit)[15]可解决这一问题。给定一个数据，FFB逐步折半搜索其最高位为1的位，从而快速检测m值。对输入数据中前部为连续0的阶，可直接调用FFB算法。对输入数据中前部为连续为1的阶，对其按位求反后转换为前部为连续0的形式，可同样调用该算法。值得注意的是，在将横纵坐标的前部均转换为连续0的形式后，仅需对其最大值进行一次置位检测即可快速获取m值。

SFO-HE算法对前m阶直接生成编码，对后(n-m)阶迭代编码，算法的执行过程如图 2所示。

图 2 SFO-HE算法的执行过程图
图选项

例1给定X=(000111)2，Y=(111101)2, 由FFB算法计算得m=3, 通过定理1计算前3阶的编码为(010101)2，第4阶的状态S4=0，然后对后3阶迭代编码，最终得编码值Z=(010101101100)2。

完整的SFO-HE算法见图 3，“~”表示按位求反。前m阶编码值可直接计算得到(第6行)，无需迭代编码。此外，本文为m由1到32分别预存储其对应的编码值，进一步提高了算法效率。

图 3 SFO-HE算法
图选项

SFO-HE算法包括置位检测和迭代编码2部分，置位检测的FFB算法的复杂度为O(lb n)，迭代编码的复杂度为O(n-m)，故算法总的复杂度为O(lb n)+ O(n-m)。复杂度随m增大而降低。

SFO-HE算法的主要空间开销在于PZMap和PSMap的开销，二者分别需要4×2×2=16 B来存储一阶编码和状态，故总的空间开销为32 B。

3 Hilbert曲线解码 3.1 状态视图的构建

解码是编码的逆过程，根据图 1，以类似编码的方式构建解码状态视图如表 3所示。

表 3 解码状态视图表
Si Z2i-1Z2i XiYi Si+1
0 00 00 1
0 01 01 0
0 10 11 0
0 11 10 3
1 00 00 0
1 01 10 1
1 10 11 1
1 11 01 2
2 00 11 3
2 01 10 2
2 10 00 2
2 11 01 1
3 00 11 2
3 01 01 3
3 10 00 3
3 11 10 0
表选项

为便于对解码算法进行描述，规定第一阶Hilbert曲线的初始状态为0(开口向下)。本文将解码状态视图设计为ZPMap和SPMap两个2维的数组。其中ZPMap存储编码值到坐标值的映射，第i阶坐标值XiYi可通过ZPMap[Si][Z2i－1Z2i]直接获取。SPMap则存储编码值到下一阶状态的映射，第(i+1)阶状态Si+1可通过SPMap[Si][Z2i－1Z2i]获取。

3.2 解码 3.2.1 SFO-HD算法

SFO-HD算法旨在跳过输入数据前部为连续的00、01、10或11的m阶，基于以下基本定理。

定理2   若前m阶编码值分别为$\overline {{Z_{2m - 1}}{Z_{2m}}} $分别为、时，设定B=2m－1，则其前m阶编码值对应的前m阶坐标为$\overline{X_{m}}=X_{1} \times B, \overline{Y_{m}}=Y_{1} \times B$，对应的第(m+1)阶的状态Sm+1分别为0(m为偶数时)或1(m为奇数时)、0、0、0(m为偶数时)或3(m为奇数时)。

证明：分4种情况分别讨论

1) $\overline {{Z_{2m - 1}}{Z_{2m}}} $为的情形。

m=1时，由Z1Z2=(00)2=0及S1=0，根据ZPMap可知$\overline {{X_1}} $=X1=0, $\overline {{Y_1}} $=Y1=0, S2=1，定理成立。

m=2时，由Z3Z4=(00)2=0, S2=1，根据ZPMap可知X2=0, Y2=0, S3=0, 故$\overline {{X_2}} $=X1X2=0, $\overline {{Y_2}} $=Y1Y2=0，定理成立。

m≥3时，第3阶输入编码值和状态与第1阶相同，故其输出坐标值和状态与第1阶相同。同理第4阶的输入和输出与第2阶相同，进入循环状态。故Sm+1=0(m为偶数)或1(m为奇数)，$\overline {{X_m}} $= ，定理成立。

2) 的情形。

m=1时，由Z1Z2=(01)2, S1=0，根据ZPMap可知$\overline {{X_1}} $=X1=0, $\overline {{Y_1}} $=Y1=1, S2=0，定理成立。

m≥2时，可见每阶的编码值和状态均与第1阶相同，则输出坐标值和下一阶状态也均与第1阶相同，进入循环状态。故可得Sm+1=0, $\overline {{X_m}} $= ，定理成立。

3) 的情形。

m=1时，由Z1Z2=(10)2, S1=0，根据ZPMap可知$\overline {{X_1}} $=X1=1, $\overline {{Y_1}} $=Y1=1, S2=0，定理成立。

m≥2时，可见每阶的输入编码值和状态均与第1阶相同，则输出坐标值和下一阶状态也均与第1阶相同，进入循环状态。故易得Sm+1=0, $\overline {{X_m}} $=Y1×B，定理成立。

4) 的情形。

m=1时，由Z1Z2=(11)2, S1=0，根据ZPMap可知$\overline {{X_1}} $=X1=1, $\overline {{Y_1}} $=Y1=0, S2=3，定理成立。

m=2时，由Z3Z4=(11)2=3, S2=3，根据ZPMap可知X2=1, Y2=0, S3=0，故$\overline {{X_2}} $=X1X2=1, $\overline {{Y_2}} $=Y1Y2=0，定理成立。

m≥3时，第3阶的编码值和状态均与第1阶相同，则输出坐标值和状态均与第1阶相同，同理第4阶的输入和输出均与第2阶相同，进入循环状态。故Sm+1=0(m为偶数)或3(m为奇数)，$\overline {{X_m}} $= ，定理成立。

故定理得证。

由定理2可见，$\overline {{X_m}} $和$\overline {{Y_m}} $分别为X1和Y1的B倍，$\overline {{X_m}} $、$\overline {{Y_m}} $和Sm+1的值可以直接计算得出，故无需迭代编码。SFO-HD算法的执行过程如图 4所示。

图 4 SFO-HD算法的执行过程图
图选项

例2给定Z=(010101101100)2, 经过FFB算法计算得m=3, 通过定理1计算前3阶的坐标$\overline {{X_3}} $=(000)2, $\overline {{Y_3}} $=(111)2，第4阶的状态S4=0，然后对后3阶迭代解码，最终得坐标X=(000111)2，Y=(111101)2。

完整的SFO-HD算法见图 5，“^”表示按位异或。对编码值开头为连续的00和11的情形，可直接对其求反后通过FFB函数检测得到m值(第3和13行)。对编码值开头为连续的01、10的情形，将其与其移位后的值进行逻辑异或，然后通过FFB函数得到m值(第6和9行)。

图 5 SFO-HD算法
图选项

考虑到解码时需针对不同的连续坐标进行相应的预处理以计算m，分支判断必不可少，故本文直接在循环中对前m阶的X和Y坐标赋值。

SFO-HD算法的主要开销在于置位检测和迭代编码2个部分，其中置位检测的FFB算法的复杂度为O(lb n)，迭代编码的复杂度为O(n-m)，故算法总的复杂度为O(lb n)+O(n-m)。随着m的增大，算法的时间复杂度逐步降低。

SFO-HD算法的主要空间开销在于ZPMap和SPMap的开销，二者分别需要4×4=16 B来存储一阶编码和状态，故总的空间开销为32 B。

4 实验 4.1 实验环境和数据集

实验硬件平台配置为：Intel i7-8750H 2.20 GHz处理器，8 GB内存。软件环境为Windows10 64位，Microsoft Visual Studio C++2010。

本文分别在模拟数据集和真实数据集下开展实验。对于模拟数据集，分别生成数据偏向左下、左上、右上和右下顶点4种情形的编码数据集，并引入偏斜阶数α来控制数据的偏斜程度，即生成的数据前部连续模式至少为α阶。为每个特定的α，各生成1百万条32阶的随机坐标数据作为实验数据。

真实数据集来自微软出租车驾驶项目T-Drive[16]，为北京10 357辆出租车的轨迹数据。原始数据为中心聚集的数据，本文选择以天安门广场为中心、经度左右各加1°、纬度上下各加0.5°范围的数据，共包含17 252 127个位置点。以天安门广场的经纬度将数据空间分割为4个区域，对每个区域做对角变换，将数据聚集于空间的4个角。

本文将生成的编码值作为解码数据集。

4.2 编码实验对比

为考察本文算法的效率，将其与领域内有代表性的文[9]、文[10]、文[12]、文[13]算法和FZF[14]算法进行对比。为便于描述，本文为上述编码算法增加后缀“-HE”，为对应的解码算法增加后缀“-HD”。

4.2.1 模拟数据集

为考察数据偏斜分布对不同算法的影响，设置分别为4、8、16、24，对这4种不同偏斜分布下的算法进行对比。图 6、7、8、9分别给出了数据向左下、左上、右上、右下偏斜时的编码时间对比。

图 6 左下偏斜编码时间
图选项

图 7 左上偏斜编码时间
图选项

图 8 右上偏斜编码时间
图选项

图 9 右下偏斜编码时间
图选项

可以看出，SFO-HE在4种偏斜分布下均具有最高的编码效率，且其编码时间均随着偏斜阶数α的增加而显著降低。左下偏斜时，SFO-HE的效率仅稍高于其他算法中效率最高的FZF-HE，这是因为FZF-HE是针对左下偏斜的情形设计的，同样可以跳过前部为0的阶，故在此情形下二者性能接近。比FZF-HE仅适合向左下偏斜的情形，SFO-HE在向4个顶点偏斜的情形下均高效，因而具有更好的偏斜分布适应性。在向左上、右上、右下偏斜时，SFO-HE的编码效率显著优于其他算法。当α=24时，SFO-HE的编码时间仅为其余算法中效率最高的文[12]-HE的37.7%、应用最广泛的文[10]-HE的8.9%，可见SFO-HE因跳过前部相同阶带来的优势。

4.2.2 真实数据集

对真实数据集T-Drive，将原始坐标分别映射为8、16、24、32阶的整数，考察它们在不同阶下的性能。考虑到文[9]-HE、文[10]-HE、文[12]-HE的效率相对较低，在此，仅将剩余3种算法进行对比，实验结果如图 10所示。

图 10 T-Drive数据集编码时间
图选项

由图 10可见，3种算法的编码时间均随阶数增加而增加，SFO-HE效率稍优于其余2种算法。尽管T-Drive数据集的平均偏斜阶数为3.4，但相比文[13]-HE，SFO-HE的编码效率仍提升8.6%。随着偏斜阶数的增加，SFO-HE的优势将进一步增大。

4.3 解码实验对比 4.3.1 模拟数据集

为考察偏斜分布的影响，设置α分别为4、8、16、24，对这4种不同偏斜分布下的算法进行比较。图 11、12、13、14分别给出了模拟数据集下数据向左下、左上、右上、右下偏斜时的解码时间对比。

图 11 左下偏斜解码时间
图选项

图 12 左上偏斜解码时间
图选项

图 13 右下偏斜解码时间
图选项

图 14 右上偏斜解码时间
图选项

可以看出，SFO-HD的解码时间均随α的增加而降低。在数据向左下偏斜时，SFO-HD与FZF-HD具有相当的编码效率，因二者均跳过前部为0的阶。除此之外，SFO-HD具有远优于其他算法的效率，具有较好的偏斜分布适应性。在数据向除左下外的3个顶点偏斜时，SFO-HD编码时间仅为其余算法中效率最高的FZF-HD的37.2%、应用最广泛的文[10]-HD的8.5%。可见SFO-HD跳过前部相同阶所带来的优势。

4.3.2 真实数据集

对T-Drive数据集，针对8、16、24、32阶的编码值考察算法在不同阶下的性能。仅对文[13]-HD、FZF-HD、SFO-HD 3种算法进行对比，结果如图 15所示。

图 15 T-Drive数据集解码时间
图选项

由图 15可见，3种算法的解码时间均随阶数增加而增加，SFO-HD效率优于其余2种算法。在数据集的平均偏斜阶数为3.4的情况下，相比其他算法中效率最高的FZF-HD，SFO-HD的解码效率可提升12.3%。

5 结论

针对二维Hilbert曲线，本文设计了高效的状态视图，并提出了一种跳过输入坐标前部为连续0或1的编码算法SFO-HE和一种跳过输入编码值前部连续00、01、10或11的解码算法SFO-HD。借助高效的位操作、置位检测算法等技术，这2个算法可避免对输入数据前部特定阶逐阶编解码。实验结果表明，本文编解码算法可提高数据向4个顶点偏斜时的编解码效率，且具有较好的数据偏斜适应性。

参考文献
[1]
SINGH H, BAWA S. A survey of traditional and mapreducebased spatial query processing approaches[J]. ACM SIGMOD Record, 2017, 46(2): 18-29. DOI:10.1145/3137586.3137590
[2]
QI J Z, LIU G L, JENSEN C S, et al. Effectively learning spatial indices[J]. Proceedings of the VLDB Endowment, 2020, 13(12): 2341-2354. DOI:10.14778/3407790.3407829
[3]
SCHRACK G, STOCCO L. Generation of spatial orders and space-filling curves[J]. IEEE Transactions on Image Processing, 2015, 24(6): 1791-1800. DOI:10.1109/TIP.2015.2409571
[4]
ZHOU L, JOHNSON C R, WEISKOPF D. Data-driven space-filling curves[J]. IEEE Transactions on Visualization and Computer Graphics, 2021, 27(2): 1591-1600. DOI:10.1109/TVCG.2020.3030473
[5]
LU J, ZHAO Y H, TAN K L, et al. Distributed density peaks cluste6ring revisited[C]//2021 IEEE 37th International Conference on Data Engineering. Chania, Greece: IEEE Press, 2021: 2352-2353.
[6]
MORTON G M. A computer oriented geodetic data base and a new technique in file sequencing[J]. Physics of Plasmas, 2015, 24(7): 159-173.
[7]
HILBERT D. Ueber die stetige abbildung einer line auf ein flächenstück[J]. Mathematische Annalen, 1970, 38(3): 459-460.
[8]
BUTZ A R. Alternative algorithm for Hilbert's space-filling curve[J]. IEEE Transactions on Computers, 1971, 100(4): 424-426.
[9]
BURKARDT J. Hilbert-curve convert between 1D and 2D coordinates of Hilbert curve[R/OL]. (2015-12-23)[2021-08-02]. http://people.math.sc.edu/Burkardt/c_src/hilbert_curve/hilbert_curve.html.
[10]
MOORE D. Fast Hilbert curve generation, sorting, and range queries[R/OL]. (2016-02-05)[2021-08-02]. https://github.com/Cheedoong/hilbert.
[11]
FISHER A J. A new algorithm for generating Hilbert curves[J]. Software: Practice and Experience, 1986, 16(1): 5-12. DOI:10.1002/spe.4380160103
[12]
李绍俊, 钟耳顺, 王少华, 等. 基于状态转移矩阵的Hilbert码快速生成算法[J]. 地球信息科学学报, 2014, 16(6): 846-851.
LI S J, ZHONG E S, WANG S H, et al. An algorithm for Hilbert ordering code based on state-transition matrix[J]. Journal of Geo-information Science, 2014, 16(6): 846-851. (in Chinese)
[13]
XCTORRES. Hilbert spatial index[R/OL]. (2017-03-30)[2021-08-02]. https://github.com/xcTorres/hilbert_spatial_index.
[14]
贾连印, 陈明鲜, 李孟娟, 等. 基于状态视图的高效Hilbert编码和解码算法[J]. 电子与信息学报, 2020, 42(6): 1494-501.
JIA L Y, CHEN M X, LI M J, et al. State view based efficient Hilbert encoding and decoding algorithms[J]. Journal of Electronics & Information Technology, 2020, 42(6): 1494-1501. (in Chinese)
[15]
ROSETTA. Find first and last set bit of a long integer[R/OL]. (2021-06-10)[2021-08-02]. http://rosettacode.org/wiki/Find_first_and_last_set_bit_of_a_long_integer.
[16]
ZHENG Y. T-Drive trajectory data sample[R/OL]. (2011-08-15)[2021-08-02]. https://www.microsoft.com/en-us/research/publication/t-drive-trajectory-data-sample/.