XZ曲线探究

最近在做时空数据库相关的项目，需要做一些空间索引的工作，于是研究了一番GeoMesa和GeoWave。发现他们专门为Polygon类型数据使用了一种叫做XZ曲线的索引。空间填充曲线对于一个GISer来说十分熟悉，但XZ曲线确实是之前闻所未闻。根据GeoMesa的介绍，找到了XZ曲线的论文。论文写得晦涩难懂，看了半天仍未完全看明白。于是又去看GeoWave中关于XZ曲线的源码，终于将XZ的基本原理大致搞懂了。

本文使用的名词解释：
cell: 传统Z曲线的网格单元
element: XZ曲线的网格单元
level: 层级
g: 最大分辨率，即最大层级
object: 对象，即需要存储的要素对象
quadrant sequences: 象限序列，即传统Z曲线的编码
interval: 一个区间(查询时使用)
interval set: 区间集合（查询时使用）

空间填充曲线

首先说一下空间填充曲线。空间填充曲线是一种将高维空间映射到一维空间的方法,经典的空间填充曲线有Z曲线、Hilbert曲线等，现在Google S2也经常被使用。Z曲线最常见的应用就是GeoHash，Hilbert是现在理论上空间邻近性最好的空间填充线曲线。关于GeoHash和Google s2可以参考高效的多维空间点索引算法 — Geohash 和 Google S2，Hilbert可以参考Hilbert曲线介绍以及代码实现。

空间填充曲线将空间进行划分为一个个大小相等的网格，一般对一个网格称为一个cell。那么，决定网格划分的数量的参数便是层级level，level越大，对空间进行的划分越精细，即level越大，cell的大小越小，数量越多。

Z曲线

Z曲线是空间填充曲线的一种。在Z曲线中，level 1将空间划分为四个cell，level 2则将level 1中的每个cell进行四等分，生成的cell的长宽是level1中cell的一半，即level 2中有16个格网。以此类推，在level n中，存在$(2^2)^n$个cell。

Z曲线更加直观地显示如图1。

图1 Z曲线

XZ曲线是根据Z曲线衍生出来的，相当于优化版的Z曲线。论文中详细阐述了为什么Z曲线需要被优化，并且XZ曲线是怎么解决Z曲线中的一些不足的。首先说明一下论文中对传统Z曲线及其编码的定义：

• Z曲线的cell按照象限序列（quadrant sequences，本文沿用论文名词）进行编码，即传统的Z曲线编码（见图2）。
• 将数据空间归一化为经纬度范围均为1的正方形，即将整个需要编码的空间作为一个边长为1的正方形，并具有x轴和y轴。则level n中每个cell的边长为$0.5^n$。

那么，作者认为，Z曲线有个致命的缺陷：

• 如果有很小的要素恰好在数据空间的正中心，则象限序列会为空（如图3）。

图2 Z曲线象限序列编码

图3 Z曲线空序列情况

XZ曲线

鉴于以上描述的Z曲线缺点，XZ曲线在Z曲线的基础上进行了优化，XZ曲线具有以下特性：

• 基本的单元（论文中称为element，本文沿用论文名词）为传统Z曲线cell向右上边长扩大一倍，临近的element之间存在50%的重叠度(如图4)。
• element的数量和cell的数量一致，即边界也向外延伸。
• XZ曲线存储要素时不作冗余存储，即单个要素只存储于一个element中。

图4 扩大的element

对于XZ曲线的一些性质和论文中给出的一些定理，本文不加证明地罗列出来，如果感兴趣可以去论文中查看证明过程。如果实在不想看这些，可以直接跳到XZ曲线的插入操作。

XZ曲线最大的一个特点就是将cell向右上扩大为element，关于element论文中有一个定义：

定义1

设一个element的左下角cell的quadrant sequences为$s$，$|s|$为其长度，则element的边长为${0.5}^{|s|-1}$。

由于一个对象（文中称为对象object，本文沿用论文名词）只存储在一个element中，则一个object最多只和一个cell的两条边相交。
则具有以下定理：

定理1

设$\lvert s\rvert$为一个quadrant sequences $s$的长度，w为一个object的宽度，h为一个object的高度。则$\lvert s\rvert$的范围为：
$l_1 \le \lvert s\rvert \lt l_2，其中 \ l_1 = \lfloor log_{0.5}(max\\{w,h\\}) \rfloor,\ l_2 = l_1 + 2$

由于不作冗余存储，一个object直接存储于一个element中，相当于使用一个element去逼近object,造成的误差会比较大。在论文中也写出了误差指标和XZ曲线的误差。
设object的实际面积为$S_{obj}$，使用一个element进行存储时element的实际面积为$S_{ele}$，则逼近误差为：
$E_{rel} = \frac{S_{ele}-S_{obj}}{S_{obj}}$
XZ曲线的最大逼近误差小于15。

设最大level为g，则对于一个level $l$的cell，其具有的level g后代数量为：
$N_{cell}(l) = 4^{g-l}$
则可求出，对于一个level $l$的cell，其具有的所有level的后代数量，是为定理2：

定理2

对于一个quadrant sequences $s$所代表的cell（level为$l$），其具有的所有level的element后代数量（包括当前对应的element，最大level为g)为：
$N_{ele} = \frac{4^{g-l+1}-1}{3}$

接下来就是最重要的一条定义，是XZ曲线的编码，称为sequence code：

定义2

对于一个quadrant sequences $s = \langle \ q_0 \ q_1 \ … \ q_i \ … \ q_{l-1} \ \rangle$，可以得到其对应的XZ曲线编码(论文中称为sequence code)为：
$C(s) = \sum_{0\le i \lt l} q_i \cdot \frac{4^{g-i}-1}{3}+1$

论文中对XZ曲线编码推论出了两个定理：

定理3

对于quadrant sequences $s$所对应的XZ编码值$C(s)$，其大小代表了$s$的字典序位置，即：
$s_1 \lt_{lex} s_2 \Longleftrightarrow C(s_1) \lt C(s_2)$

定理4

从quadrant sequences映射到数值型编码的方法中，sequence code具有最好的临近性。即相近的两个quadrant sequences对应的sequence code值最相近。

XZ曲线的插入操作

XZ曲线的插入操作实际上就是获取一个object对应的XZ曲线编码。在GIS领域里，一般为了方便对不规则形状进行管理，都会将object加上一个BoundingRectangle，可以想象成一个能够将这个object装下的最小矩形。这个矩形存储object的详细信息，需要做拓扑计算等操作的时候可以使用BoundingRectangle进行一个快速的预计算。
首先，需要确定object应该所处的层级，以确保一个element中能装下这个object的BoundingRectangle。
根据定理1，设：
$l_1 = \lfloor log_{0.5}(max\\{w,h\\})\rfloor$
其中，$w$为BoundingRectangle的宽度，$h$为BoundingRectangle的高度，一个object所在的层级为$l_1$或$l_1+1$。
设$x_l$为BoundingRectangle的横轴（x轴）的下界，$y_l$为object的纵轴（y轴）的下界。设$w_2 = {0.5}^{l_1+1}$。则可根据以下两个式子可判断object的层级$l$：
$x_l + w \le \lfloor \frac{x_l}{w_2} \rfloor \cdot w_2 + 2w_2$
$y_l + h \le \lfloor \frac{y_l}{w_2} \rfloor \cdot w_2 + 2w_2$
注：这个不等式和论文上不一致，经过讨论发现应该是论文出错，此处公式使用的是GeoWave源码中的公式
若以上两式同时成立，则object的层级$l = l_1 + 1$，否则$l = l_1$

得到层级$l$后，通过对空间的不断四分，得到BoundingRectangle左下角点所在的第$l$层的Z曲线编码(quadrant sequences)，然后，根据定义2，就可以得到XZ曲线编码(sequence code)。

XZ曲线的查询操作

XZ曲线查询得到的结果是一个区间集合（interval set，本文沿用论文名词），interval set中的一个interval代表一个element中所有后代element的集合，即interval表示的是sequence code范围。

关于XZ曲线的查询生成interval set，论文中表述得较少并且不明确，GeoWave源码也没有完全按照论文中的方法来实现。本着实事求是的原则，以下按照GeoWave源码的方法讲述XZ曲线查询操作。

查询的方法为：给定一个查询框，如果查询框和包含某个element，则将该element对应的interval（即其所有后代element）加入到interval set中；如果查询框与某个element相交，则将该element作为一个interval（这个interval中只有这个element对象）加入到interval set中，并将这个element分裂为其子element，再和查询框进行拓扑关系判断，直至interval set中的数量达到给定的阈值或分裂至最大层级。得到interval set后，将其中有交集的interval进行合并，生成最终的interval set。

查询的伪代码如下：

Queue<Element> remaining; // 待处理队列
List<Interval> intervalSet; // 最终结果
remaining.add(LevelOneElement); // 加入第一层的element
remaining.add(LevelTerminator); // 加入一层的终止条件

while (level < g && !remaining.empty() && intervalSet.size() < maxSize) {
    Element next = remaining.poll();
    if (next == LevelTerminator) {
        if (!remaining.isEmpty()) {
            level = level + 1;
            remaining.add(LevelTerminator);
        }
    }
    else {
        if (query.contain(next)) {
            intervalSet.add(interval(next));
        }
        else if (query.intersect(next)) {
            intervalSet.add(next);
            remaining.add(next.children());
        }
    }
}

// 兜底策略，将未处理完的element处理完
while (!remaining.empty()) {
    Element next = remaining.poll();
    if (next == LevelTerminator) {
        level = level + 1;
    }
    else {
        intervalSet.add(interval(next));
    }
}

mergeIntervalsWhichAreIntersected(intervalSet);

总结

XZ曲线的论文在1999年就发表了，是一篇会议论文，但是基本上没什么名气，用的人也很少，目前来看除了GeoMesa/GeoWave用过外，就没有什么知名开源项目使用过了，所以现在能获取的相关信息很少（感觉有一个原因是论文写得太晦涩难懂了），只能通过原始的论文和GeoWave的代码来一步步搞清楚原理。
论文作者在实验章节基于ORACLE-8做了一个实验（然而这个代码没公开），测试XZ曲线和传统Z曲线查询时磁盘页的访问数量，结果表明XZ曲线的磁盘页访问数量远小于Z曲线。然而我感觉虽然Z曲线会冗余存储，但是在查询速度上可能是要胜于XZ曲线的。
最后，在好不容易看完这篇论文后，我们又调研了其他的分布式时空数据库，然而没有发现其他使用XZ曲线的，所以我也没有测试XZ曲线的性能。如果看到了这篇博客并且有兴趣的大佬，可以尝试把GeoWave的代码clone下来后跑一跑他们XZ曲线相关的代码（项目地址在博客第一段已经给出）。