针对三维路网轨迹数据（经度、维度、时间戳），需要将这些三维数据转化为有序的一维数据，实现空间数据的集中存储，以减少查询IO开销。常见的空间填充曲线包括Z曲线、希尔伯特曲线和XZ-Ordering。空间填充曲线具有出色的邻近保持属性，使其成为线性化数据对象的理想选择。本文将涉及Z曲线和希尔伯特曲线，并着重讲解下希尔伯特曲线。

空间填充曲线可以分为两类：单调空间填充曲线和非单调空间填充曲线。两者核心区别在于空间上的几何顺序是否能直接对应到一维映射后的数值顺序。如果一个空间填充曲线的映射函数满足：只要空间上 ，那么映射后的一维数值 恒成立，则该曲线是单调的。单调空间填充曲线，如Z曲线可以快速定位搜索范围，无论查询窗口有多大，只需要计算两个顶点的映射值，就能确定搜索范围，计算复杂度低。非单调空间填充曲线，如希尔伯特曲线，拥有更好的邻近性保持，空间局部性保持能力强，在一维序列上相邻的两个点，在原始多维空间中总是非常接近的。

Z曲线（Z-order Curve）是一种将多维空间数据映射到一维空间的方法。把地图上的二维坐标映射为一个单独的整数索引。在遍历网格节点时，它的轨迹呈现出字母 “Z” 的形状，被称为Z曲线。如下图所示，Z曲线递归地将空间分成四个子空间，直到达到最大递归次数r，其中四个子空间分别按照图所示方式从0到3编号。最大分辨率控制着最小网格的大小。如果没有达到最大分辨率，则继续划分每一个子空间，并依次递归编码。

希尔伯特曲线是一种能填充满一个平面正方形的分形曲线，该曲线把一个正方形空间不断划分为 4 个子空间，再把小正方形的中心点连接起来。在希尔伯特曲线的编码映射中，使用 U 字型来访问每个空间，对分成的 4 个子空间也同样使用 U 字形访问，但要调整 U 字形的朝向使得相邻的空间能够衔接起来。希尔伯特曲线示例如下图所示。

相比于其他空间填充曲线（如Z-order曲线），希尔伯特曲线的主要优点是空间特征保存效果好，多维空间中的连续区域映射到一维希尔伯特值后，通常会分布在连续邻近的区间内，邻近保持性强。

下面我将介绍基于希尔伯特曲线的排序方法。这是将高维数据映射为一维序号的重要步骤。

定义1（方向模式）给定一个有四个单元格的子网格SG ，它的方向模式被设为一个2*2的矩阵dp[i,j]∈{0,1,2,3}(0<=i,j<=1) ，表示子网格内的四个单元格被希尔伯特曲线访问的局部先后顺序。

定义2（父方向模式）给定一个子网格SG和其方向模式dp，记SG内部的一个单元格gc的父方向模式为dp，记作pp(gc)。

图1 一个四单元格网格及其序列号

图1展示了一个包含四个单元格的网格，其初始序列号为{0,1,2,3}，这与一阶希尔伯特曲线一致。网

格的方向模式为

首先，将所有的轨迹样本放置在一个网格中，该网格的空间填充希尔伯特曲线阶数为零。然后需要将这个网格划分为4个子网格，希尔伯特曲线阶数为一。每次进行网格的划分，希尔伯特曲线阶段加一。假设一个由4个单元格组成的子网格被一条阶数为z的希尔伯特曲线穿过，其方向模式为dp(SG)，如果每一个单元格gc∈SG被进一步划分为一个四单元格的网格，希尔伯特曲线阶段为z+1，这4个单元格内部的曲线方向是由它在父格中的位置决定的。

局部希尔伯特序号为1或者2的单元格对应的子网格的方向模式与dp(SG)相同；如果位于父网格遍历顺序的起始位置（序号为0）或终止位置（序号为3），子网格的方向模式需要进行变换，通过沿着“不经过当前起始（或结束）点的对角线”交换数值来获得子方向模式。因此，采用算法1来实现上述逻辑。

算法1：子方向模式判断numPtn
输入：父方向模式pp，当前单元格在父模式中的序号nm
输出：子网格的方向模式
1.if nm=1 or nm=2 then
2.   return pp;
3.else if nm=0 then
4.  [i1,j1]←pp局部序号为1的单元格坐标;
5.  [i3,j3]←pp局部序号为3的单元格坐标;
6.  swap pp[i1,j1]  and pp[i3,j3];
7.  return pp;
8.else if nm=3 then
9.   [i0,j0]←pp局部序号为0的单元格坐标;
10.  [i2,j2]←pp局部序号为2的单元格坐标;
11.  swap pp[i0,j0]  and pp[i2,j2];
12.  return pp;
13.end if

图2 局部序号划分网格

图2展示了图1的细化网格，希尔伯特曲线阶数为2，已知图1的方向模式为，当进行网格的细化划分时，根据父方向模式pp=，根据算法1描述，序号1和2的单元格进行划分，不改变方向模式，所以左上和右上单元格的方向模式仍为，而序号0和3对应的单元格需要交换对角线值以获得新方向模式，因此左下单元格的方向模式为，右下单元格的方向模式为。图3展示了二阶希尔伯特曲线在细化网格中的编号。

图3 最终希尔伯曲线序号图

对原始样本数据进行排序，以此确定每个样本的序列号和希尔伯特曲线的最大阶数。当单元格内样本数量大于样本阈值μ时，需要对该单元格进行划分，完整的划分过程如算法2所示，主要思路为当初始的轨迹数据集包含μ个样本数据，则直接根据输入的起始序列号来确定样本数据的序列号，如果初始的轨迹数据集包含更多样本，则根据初始方向模式将时间-空间二维区域划分为4个单元格。左下、左上、右上和右下单元格分别对应D_l^b，D_l^t，D_r^t，D_r^b ，每个单元格分别计算其方向模式和细化子网格中的起始序列号，并递归调用算法进行进一步划分。

算法2：网格划分divideConquerSort
输入：空间维度范围[p1,p2]，时间维度范围[t1,t2]，轨迹数据集D，初始方向模式idp，起始序列号ssn，样本阈值μ
输出：轨迹数据集中每个样本的序列号和最大希尔伯特曲线阶数
1.if |D|≤μ then
2.      if |D|>0 then
3.       for each sample p in D do
4.           p.sn←ssn;
5.       end for
6.      ssn←ssn+1;
7.   end if
8.   return 0;
9.else if |D|>μ then
10.   基于[p1,p2]的中点pm和[t1,t2]的中点tm将D划分为{Dlb，Dlt，Drt，Drb};
11.   for  k←0; k<4;k++ do
12.           i,j←arg⁡idp[i,j]=k;  //确定当前k对应哪个物理象限(i, j)
13.           targetD←Dij;    //获取该象限的数据集
14.           newidp←numPtn(idp,k);  //获取该象限的方向模式
15.          if |targetD|>0 then //只有当前区域有数据时，进行递归划分
16.           depth ←divideConquerSort(targetD,newidp,ssn);
17.           max⁡_depth ←max⁡(depth,max⁡_depth);
18.       end if
19.   end for
20.   return 1+max⁡_depth;
21.end if

算法2是一个典型的分治法算法，结构上类似于四叉树的构造过程，对空间进行递归划分并赋予希尔伯特序号，如果数据在空间分布比较均匀，每次划分都将数据N大致四等分，树的高度为log_4^N，可得时间复杂度为O(NlogN)。

图4 网格划分实例

例：根据算法2，设定样本阈值μ=2，起始序列号ssn=0 ，初始方向模式idp=，如图4所示，时间-空间平面内一共有10个时空样本数据，现在进行网格划分与数据排序。首先，初始平面内尚未进行划分，10个样本数量超过了设定的样本阈值μ=2，根据算法10行，通过中点将平面划分为四个区域，构成一个有四个单元格的网格。接着根据算法13、15行分别统计四个区域的时空样本数量并与样本阈值μ比较，左下、左上、右上、右下区域分别包含3、0、4、3个样本，因此除了左上区域之外，其他三个区域都需要进一步递归划分。算法14行调用了算法1numPtn计算子区域的方向模式，以此根据方向模式计算子区域内的希尔伯特序号。对于左下区域，由于样本数量仍然超过样本阈值μ，所以还需要进一步划分直到每个单元格内样本数量小于μ。根据序列号ssn=0和方向模式，按顺序对单元格内的样本点赋值，样本点a、b处于同一单元格内被赋值为0，下一单元格内的样本点c被赋值为1。按照上述操作，对每个单元格内的样本赋值并排序。由此得到最终的序号排序{[a,b],[c],[d],[e,f],[g],[h,i],[j]}={0,1,2,3,4,5,6}。