【转载】高效的多维空间点索引算法 — Geohash 和 Goo
转载一篇文章，因为它解决了一直以来的疑惑，之前做过数据存储的项目，虽然采用了曲线方式，但却不理解其究竟，看了它之后懂了。原文地址
高效的多维空间点索引算法 — Geohash 和 Google S2
每天我们晚上加班回家，可能都会用到滴滴或者共享单车。打开 app 会看到如下的界面：
app 界面上会显示出自己附近一个范围内可用的出租车或者共享单车。假设地图上会显示以自己为圆心，5公里为半径，这个范围内的车。如何实现呢？最直观的想法就是去数据库里面查表，计算并查询车距离用户小于等于5公里的，筛选出来，把数据返回给客户端。
这种做法比较笨，一般也不会这么做。为什么呢？因为这种做法需要对整个表里面的每一项都计算一次相对距离。太耗时了。既然数据量太大，我们就需要分而治之。那么就会想到把地图分块。这样即使每一块里面的每条数据都计算一次相对距离，也比之前全表都计算一次要快很多。
我们也都知道，现在用的比较多的数据库 MySQL、PostgreSQL 都原生支持 B+ 树。这种数据结构能高效的查询。地图分块的过程其实就是一种添加索引的过程，如果能想到一个办法，把地图上的点添加一个合适的索引，并且能够排序，那么就可以利用类似二分查找的方法进行快速查询。
问题就来了，地图上的点是二维的，有经度和纬度，这如何索引呢？如果只针对其中的一个维度，经度或者纬度进行搜索，那搜出来一遍以后还要进行二次搜索。那要是更高维度呢？三维。可能有人会说可以设置维度的优先级，比如拼接一个联合键，那在三维空间中，x，y，z 谁的优先级高呢？设置优先级好像并不是很合理。
本篇文章就来介绍2种比较通用的空间点索引算法。
Geohash 是一种地理编码，由 Gustavo Niemeyer 发明的。它是一种分级的数据结构，把空间划分为网格。Geohash 属于空间填充曲线中的 Z 阶曲线（Z-order curve）的实际应用。
何为 Z 阶曲线？
上图就是 Z 阶曲线。这个曲线比较简单，生成它也比较容易，只需要把每个 Z 首尾相连即可。
Z 阶曲线同样可以扩展到三维空间。只要 Z 形状足够小并且足够密，也能填满整个三维空间。
说到这里可能读者依旧一头雾水，不知道 Geohash 和 Z 曲线究竟有啥关系？其实 Geohash算法 的理论基础就是基于 Z 曲线的生成原理。继续说回 Geohash。
Geohash 能够提供任意精度的分段级别。一般分级从 1-12 级。
字符串长度cell 宽度cell 高度
1≤5,000km×5,000km
2≤1,250km×625km
3≤156km×156km
4≤39.1km×19.5km
5≤4.89km×4.89km
6≤1.22km×0.61km
7≤153m×153m
8≤38.2m×19.1m
9≤4.77m×4.77m
10≤1.19m×0.596m
11≤149mm×149mm
12≤37.2mm×18.6mm
还记得引语里面提到的问题么？这里我们就可以用 Geohash 来解决这个问题。
我们可以利用 Geohash 的字符串长短来决定要划分区域的大小。这个对应关系可以参考上面表格里面 cell 的宽和高。一旦选定 cell 的宽和高，那么 Geohash 字符串的长度就确定下来了。这样我们就把地图分成了一个个的矩形区域了。
地图上虽然把区域划分好了，但是还有一个问题没有解决，那就是如何快速的查找一个点附近邻近的点和区域呢？
Geohash 有一个和 Z 阶曲线相关的性质，那就是一个点附近的地方(但不绝对) hash 字符串总是有公共前缀，并且公共前缀的长度越长，这两个点距离越近。
由于这个特性，Geohash 就常常被用来作为唯一标识符。用在数据库里面可用 Geohash 来表示一个点。Geohash 这个公共前缀的特性就可以用来快速的进行邻近点的搜索。越接近的点通常和目标点的 Geohash 字符串公共前缀越长（但是这不一定，也有特殊情况，下面举例会说明）
Geohash 也有几种编码形式，常见的有2种，base 32 和 base 36。
Decimal0123456789101112131415
Base 320123456789bcdefg
Decimal16171819202122232425262728293031
Base 32hjkmnpqrstuvwxyz
base 36 的版本对大小写敏感，用了36个字符，“23456789bBCdDFgGhHjJKlLMnNPqQrRtTVWX”。
Decimal0123456789101112131415161718
Base 3623456789bBCdDFgGhHj
Decimal1920212223242526272829303132333435
Base 36JKILMnNPqQrRtTVWX
接下来的举例以 base-32 为例。举个例子。
上图是一个地图，地图中间有一个美罗城，假设需要查询距离美罗城最近的餐馆，该如何查询？
第一步我们需要把地图网格化，利用 geohash。通过查表，我们选取字符串长度为6的矩形来网格化这张地图。
经过查询，美罗城的经纬度是[31.1932993, 121.43960190000007]。
先处理纬度。地球的纬度区间是[-90,90]。把这个区间分为2部分，即[-90,0)，[0,90]。31.1932993位于(0,90]区间，即右区间，标记为1。然后继续把(0,90]区间二分，分为[0,45)，[45,90]，31.1932993位于[0,45)区间，即左区间，标记为0。一直划分下去。
左区间中值右区间二进制结果
-900901
045900
022.5451
22.533.75450
22.528.12533.751
28.12530.937533.751
30.937532.3437533.750
30.937531.64062532.343750
30.937531.289062531.6406250
30.937531.113281231.28906251
31.113281231.201171831.28906250
31.113281231.157226531.20117181
31.157226531.179199231.20117181
31.179199231.190185531.20117181
31.190185531.195678631.20117180
再处理经度，一样的处理方式。地球经度区间是[-180,180]
左区间中值右区间二进制结果
-18001801
0901801
901351800
90112.51351
112.5123.751350
112.5118.125123.751
118.125120.9375123.751
120.9375122.34375123.750
120.9375121.640625122.343750
120.9375121.289062121.6406251
121.289062121.464844121.6406250
121.289062121.376953121.4648441
121.376953121.420898121.4648441
121.420898121.442871121.4648440
121.420898121.431885121.4428711
纬度产生的二进制是101011000101110，经度产生的二进制是110101100101101，按照“偶数位放经度，奇数位放纬度”的规则，重新组合经度和纬度的二进制串，生成新的：111001100111100000110011110110，最后一步就是把这个最终的字符串转换成字符，对应需要查找 base-32 的表。11100 11001 11100 00011 00111 10110转换成十进制是 28 25 28 3 7 22，查表编码得到最终结果，wtw37q。
我们还可以把这个网格周围8个各自都计算出来。
从地图上可以看出，这邻近的9个格子，前缀都完全一致。都是wtw37。
如果我们把字符串再增加一位，会有什么样的结果呢？Geohash 增加到7位。
当Geohash 增加到7位的时候，网格更小了，美罗城的 Geohash 变成了 wtw37qt。
看到这里，读者应该已经清楚了 Geohash 的算法原理了。咱们把6位和7位都组合到一张图上面来看。
可以看到中间大格子的 Geohash 的值是 wtw37q，那么它里面的所有小格子前缀都是 wtw37q。可以想象，当 Geohash 字符串长度为5的时候，Geohash 肯定就为 wtw37 了。
接下来解释之前说的 Geohash 和 Z 阶曲线的关系。回顾最后一步合并经纬度字符串的规则，“偶数位放经度，奇数位放纬度”。读者一定有点好奇，这个规则哪里来的？凭空瞎想的？其实并不是，这个规则就是 Z 阶曲线。看下图：
x 轴就是纬度，y轴就是经度。经度放偶数位，纬度放奇数位就是这样而来的。
最后有一个精度的问题，下面的表格数据一部分来自 Wikipedia。
Geohash 字符串长度纬度经度纬度误差经度误差km误差
123±23±23±2500
255±2.8±5.6±630
378±0.70±0.70±78
41010±0.087±0.18±20
51213±0.022±0.022±2.4
61515±0.0027±0.0055±0.61
71718±0.00068±0.00068±0.076
82020±0.000085±0.00017±0.019
92223
102525
112728
123030
到此，读者应该对 Geohash 的算法都很明了了。接下来用 Go 实现一下 Geohash 算法。
packagegeohashimport("bytes")const(BASE32="0123456789bcdefghjkmnpqrstuvwxyz"MAX_LATITUDEfloat64=90MIN_LATITUDEfloat64= -90MAX_LONGITUDEfloat64=180MIN_LONGITUDEfloat64= -180)var(bits = []int{16,8,4,2,1}base32 = []byte(BASE32))typeBoxstruct{MinLat,MaxLatfloat64//纬度MinLng,MaxLngfloat64//经度}func(this*Box)Width()float64{returnthis.MaxLng- this.MinLng}func(this*Box)Height()float64{returnthis.MaxLat- this.MinLat}//输入值：纬度，经度，精度(geohash的长度)//返回geohash, 以及该点所在的区域funcEncode(latitude,longitudefloat64,precisionint) (string, *Box) {vargeohashbytes.BuffervarminLat,maxLatfloat64= MIN_LATITUDE, MAX_LATITUDEvarminLng,maxLngfloat64= MIN_LONGITUDE, MAX_LONGITUDEvarmidfloat64=0bit,ch,length,isEven:=0,0,0,trueforlength < precision {ifisEven {ifmid = (minLng + maxLng) /2; mid < longitude {ch |= bits[bit]minLng = mid}else{maxLng = mid}}else{ifmid = (minLat + maxLat) /2; mid < latitude {ch |= bits[bit]minLat = mid}else{maxLat = mid}}isEven = !isEvenifbit <4{bit++}else{geohash.WriteByte(base32[ch])length, bit, ch = length+1,0,0}}b:=&Box{MinLat: minLat,MaxLat: maxLat,MinLng: minLng,MaxLng: maxLng,}returngeohash.String(), b}
Geohash 的优点很明显，它利用 Z 阶曲线进行编码。而 Z 阶曲线可以将二维或者多维空间里的所有点都转换成一维曲线。在数学上成为分形维。并且 Z 阶曲线还具有局部保序性。
Z 阶曲线通过交织点的坐标值的二进制表示来简单地计算多维度中的点的z值。一旦将数据被加到该排序中，任何一维数据结构，例如二叉搜索树，B树，跳跃表或（具有低有效位被截断）哈希表 都可以用来处理数据。通过 Z 阶曲线所得到的顺序可以等同地被描述为从四叉树的深度优先遍历得到的顺序。
这也是 Geohash 的另外一个优点，搜索查找邻近点比较快。
Geohash 的缺点之一也来自 Z 阶曲线。
Z 阶曲线有一个比较严重的问题，虽然有局部保序性，但是它也有突变性。在每个 Z 字母的拐角，都有可能出现顺序的突变。
看上图中标注出来的蓝色的点点。每两个点虽然是相邻的，但是距离相隔很远。看右下角的图，两个数值邻近红色的点两者距离几乎达到了整个正方形的边长。两个数值邻近绿色的点也达到了正方形的一半的长度。
Geohash 的另外一个缺点是，如果选择不好合适的网格大小，判断邻近点可能会比较麻烦。
看上图，如果选择 Geohash 字符串为6的话，就是蓝色的大格子。红星是美罗城，紫色的圆点是搜索出来的目标点。如果用 Geohash 算法查询的话，距离比较近的可能是 wtw37p，wtw37r，wtw37w，wtw37m。但是其实距离最近的点就在 wtw37q。如果选择这么大的网格，就需要再查找周围的8个格子。
如果选择 Geohash 字符串为7的话，那变成黄色的小格子。这样距离红星星最近的点就只有一个了。就是 wtw37qw。
如果网格大小，精度选择的不好，那么查询最近点还需要再次查询周围8个点。
在介绍第二种多维空间点索引算法之前，要先谈谈空间填充曲线(Space-filling curve)和分形。
解决多维空间点索引需要解决2个问题，第一，如何把多维降为低维或者一维？第二，一维的曲线如何分形？
在数学分析中，有这样一个难题：能否用一条无限长的线，穿过任意维度空间里面的所有点？
在1890年，Giuseppe Peano 发现了一条连续曲线，现在称为 Peano 曲线，它可以穿过单位正方形上的每个点。他的目的是构建一个可以从单位区间到单位正方形的连续映射。 Peano 受到 Georg Cantor 早期违反直觉的研究结果的启发，即单位区间中无限数量的点与任何有限维度流型（manifold）中无限数量的点，基数相同。 Peano 解决的问题实质就是，是否存在这样一个连续的映射，一条能填充满平面的曲线。上图就是他找到的一条曲线。
一般来说，一维的东西是不可能填满2维的方格的。但是皮亚诺曲线恰恰给出了反例。皮亚诺曲线是一条连续的但处处不可导的曲线。
皮亚诺曲线的构造方法如下：取一个正方形并且把它分出九个相等的小正方形，然后从左下角的正方形开始至右上角的正方形结束，依次把小正方形的中心用线段连接起来；下一步把每个小正方形分成九个相等的正方形，然后上述方式把其中中心连接起来……将这种操作手续无限进行下去，最终得到的极限情况的曲线就被称作皮亚诺曲线。
皮亚诺对区间[0，1]上的点和正方形上的点的映射作了详细的数学描述。实际上，正方形的这些点对于，可找到两个连续函数 x = f(t) 和 y = g(t)，使得 x 和 y 取属于单位正方形的每一个值。
一年后，即1891年，希尔伯特就作出了这条曲线，叫希尔伯特曲线（Hilbert curve）。
上图就是1-6阶的希尔伯特曲线。具体构造方式在下一章再说。
上图是希尔伯特曲线填充满3维空间。
之后还有很多变种的空间填充曲线，龙曲线(Dragon curve)、 高斯帕曲线(Gosper curve)、Koch曲线(Koch curve)、摩尔定律曲线(Moore curve)、谢尔宾斯基曲线(Sierpiński curve)、奥斯古德曲线(Osgood curve)。这些曲线和本文无关，就不详细介绍了。
在数学分析中，空间填充曲线是一个参数化的注入函数，它将单位区间映射到单位正方形，立方体，更广义的，n维超立方体等中的连续曲线，随着参数的增加，它可以任意接近单位立方体中的给定点。除了数学重要性之外，空间填充曲线也可用于降维，数学规划，稀疏多维数据库索引，电子学和生物学。空间填充曲线的现在被用在互联网地图中。
皮亚诺曲线的出现，说明了人们对维数的认识是有缺陷的，有必要重新考察维数的定义。这就是分形几何考虑的问题。在分形几何中，维数可以是分数叫做分维。
多维空间降维以后，如何分形，也是一个问题。分形的方式有很多种，这里有一个列表，可以查看如何分形，以及每个分形的分形维数，即豪斯多夫分形维(Hausdorff fractals dimension)和拓扑维数。这里就不细说分形的问题了，感兴趣的可以仔细阅读链接里面的内容。
接下来继续来说多维空间点索引算法，下面一个算法的理论基础来自希尔伯特曲线，先来仔细说说希尔伯特曲线。
希尔伯特曲线一种能填充满一个平面正方形的分形曲线（空间填充曲线），由大卫·希尔伯特在1891年提出。
由于它能填满平面，它的豪斯多夫维是2。取它填充的正方形的边长为1，第n步的希尔伯特曲线的长度是2^n - 2^(-n)。
一阶的希尔伯特曲线，生成方法就是把正方形四等分，从其中一个子正方形的中心开始，依次穿线，穿过其余3个正方形的中心。
二阶的希尔伯特曲线，生成方法就是把之前每个子正方形继续四等分，每4个小的正方形先生成一阶希尔伯特曲线。然后把4个一阶的希尔伯特曲线首尾相连。
三阶的希尔伯特曲线，生成方法就是与二阶类似，先生成二阶希尔伯特曲线。然后把4个二阶的希尔伯特曲线首尾相连。
n阶的希尔伯特曲线的生成方法也是递归的，先生成n-1阶的希尔伯特曲线，然后把4个n-1阶的希尔伯特曲线首尾相连。
看到这里可能就有读者有疑问了，这么多空间填充曲线，为何要选希尔伯特曲线？
因为希尔伯特曲线有非常好的特性。
首先，作为空间填充曲线，希尔伯特曲线可以对多维空间有效的降维。\
上图就是希尔伯特曲线在填满一个平面以后，把平面上的点都展开成一维的线了。可能有人会有疑问，上图里面的希尔伯特曲线只穿了16个点，怎么能代表一个平面呢？
当然，当n趋近于无穷大的时候，n阶希尔伯特曲线就可以近似填满整个平面了。
当n阶希尔伯特曲线，n趋于无穷大的时候，曲线上的点的位置基本上趋于稳定。举个例子：
上图左边是希尔伯特曲线，右边是蛇形的曲线。当n趋于无穷大的时候，两者理论上都可以填满平面。但是为何希尔伯特曲线更加优秀呢？
在蛇形曲线上给定一个点，当n趋于无穷大的过程中，这个点在蛇形曲线上的位置是时刻变化的。
这就造成了点的相对位置始终不定。再看看希尔伯特曲线，同样是一个点，在n趋于无穷大的情况下：
从上图可以看到，点的位置几乎没有怎么变化。所以希尔伯特曲线更加优秀。
希尔伯特曲线是连续的，所以能保证一定可以填满空间。连续性是需要数学证明的。具体证明方法这里就不细说了，感兴趣的可以点文章末尾一篇关于希尔伯特曲线的论文，那里有连续性的证明。
接下来要介绍的谷歌的 S2 算法就是基于希尔伯特曲线的。现在读者应该明白选择希尔伯特曲线的原因了吧。
