HBase^[3]是开源的面向列族的分布式NoSQL数据库，构建在分布式文件系统HDFS(Hadoop distributed file system)之上，实现了基于Key-Value数据模型的数据存储方式。HBase中的数据被组织成表的形式，表中的每一行包含一个主键(Rowkey)和多个列族，每个列族又包含多个列。HBase为每个表的主键创建了分布式B+树索引，给定Rowkey可以快速访问对应的数据行。HBase还支持针对Rowkey的Scan查询，即指定Rowkey的起始值和结束值，返回Rowkey在该范围内的所有数据行。由于没有建立索引，对HBase表中非主键的属性列进行查询时，必须进行全表扫描，所以检索效率很低。HBase支持数据的高速插入，随着数据量的增加，HBase表将会自动分裂成多个子表(Region)，不同的Region一般由不同的服务器存储和管理。

HBase仅提供专用的数据访问接口，不支持标准SQL接口。为此，Apache的开源项目Phoenix^[10]在HBase平台上实现了基于标准SQL语言的数据访问接口支持。Phoenix将应用层发出的SQL操作请求转换为HBase的本地应用程序接口(application programming interface, API)调用，并利用HBase协处理器技术尽可能将查询处理放在HBase集群上并行执行。Phoenix还为HBase增加了二次索引、用户自定义函数、服务器端数据聚集和排序、数据视图、查询结果分页、多表连接查询和跨表事务支持等多种现代关系型数据库具有的高级数据管理机制。

存储二维空间数据的一种常见的方法是利用空间填充曲线将高维数据一维化，再用传统的一维方式建立索引。空间填充曲线能够最大限度地使高维空间上临近的两个点在一维空间上也保持临近。空间数据一维化后可以获得一个唯一的编码值，称为Geohash编码。常用的二维空间填充曲线包括Z曲线^[11]、Hilbert曲线^[12]和Peano曲线等。其中，Z曲线虽然在保持空间临近性上略逊于Hilbert曲线，但由于更易于实现编码和解析，所以应用更为广泛。本文也采用Z曲线实现Geohash编码。

不失一般性，假定要编码的整个二维空间是一个正方形。在该空间上基于Z曲线实现Geohash编码的方法如下：首先根据编码阶数l将整个二维空间的x轴和y轴跨度平分为2$^{l}$份，整个空间被均匀划分为2$^{2l}$个单元格。每个单元格的x轴和y轴坐标都可以用一个l位的二进制串表示。将任意单元格的横轴和纵轴坐标的二进制串按位交叉叠加在一起，形成一个2l位的二进制串，就是该单元格在一维空间的编码，也是所有落在该单元格内的空间点对象的Geohash编码。图 1给出了阶数为1、2、3的Geohash编码方法。

图  1  Geohash编码过程

Figure 1.  Geohash Coding Procedure

GS-Phoenix是面向地理空间数据的存储管理系统，其主要设计目的是为HBase增加带空间扩展功能的SQL语言访问能力，使系统支持地理空间点类型的数据定义，支持关于矩形、多边形和圆形的空间范围查询、kNN查询，支持Distance、Buffer等常用的空间函数计算。

GS-Phoenix的系统构架设计如图 2所示，主要由客户端和服务器端两个部分组成。应用层的SQL访问请求首先被客户端的SQL解析器截获。如果是不包含空间查询的普通SQL语句，则直接将该请求转发给Phoenix的客户端API；否则将二维空间查询转化为多个一维空间的区间查询，并生成Phoenix支持的SQL访问。GS-Phoenix将改写后的SQL语句转发给Phoenix的客户端API，获得查询结果后，经过整理再返回给应用层。GS-Phoenix客户端还将利用数据分布统计模块获得的数据空间分布统计信息优化查询方案。在服务器端，数据被存储到HBase表中，并由空间属性字段创建了主键索引或二次索引。客户端的空间查询请求被转化为针对HBase表的多个Scan操作。Scan操作粗选的记录集合再利用空间UDF函数进行精细过滤，因此空间查询中的主要运算任务被放置在服务器端，并且避免了不必要的数据传输。本文主要利用开源的GeoTools工具库^[13]实现UDF过滤中所需的空间关系运算。

图  2  GS-Phoenix系统构架

Figure 2.  Architecture of GS-Phoenix

GS-Phoenix根据用户的数据模型定义，利用Phoenix支持的数据定义语句(data definition language, DDL) 创建对应的HBase数据表，数据模型中的各个字段都有HBase表中的列与之对应。空间位置字段的经度和纬度坐标分别以浮点数类型存储到HBase表的两个不同列中。GS-Phoenix将自动为HBase表增加一个名为"geohash"的属性列，用于存储空间位置的Geohash编码。如果用户指定空间位置字段作为数据表的主键，则Geohash列被作为HBase表的行键或者与其他列合并后作为行键。如果用户指定其他属性字段作为主键，则利用Phoenix的CREATE INDEX命令创建关于Geohash列的二次索引。因此，GS-Phoenix可以同时支持对空间属性和非空间属性的有效查询。

假定数据表mytable包含entity_id、lo、la、c1、c1、c3共6列，其中entity_id为对象标识，lo和la字段分别存储了每条记录所绑定空间点的经度和纬度坐标，c1、c2、c3是其他非空间属性字段。GS-Phoenix调用如下DDL语句创建表mytable，并将(lo, la)对应的Geohash嵌入到表的主键：

CREATE TABLE mytable (geohash CHAR(2) primary key, entity_id CHAR(10), lo DOUBLE, la DOUBLE, c1 VARCHAR, c2 VARCHAR, c3, CONSTRAINT pk PRIMARY KEY (geohash, entity_id))

GS-Phoenix实现空间范围查询的基本思路是：先将矩形区域转化为Geohash一维空间的多个连续区间，然后在一维空间实施区间查询。因为在数据表或二次索引表中已经关于Geohash列建立了分布式B+树索引，所以针对Geohash列的区间查询具有很高的效率。

假定要针对§2.2定义的表mytable作空间查询，其中Geohash列采用图 3所示的6位Geohash编码。图 3中，整个二维空间被划分为64个单元格，每个单元格对应的二进制和十进制形式的Geohash编码标注在单元格内部。红色粗线所示的二维矩形为查询区域，则该区域可以转化为[16, 31]、[48, 49]和[52, 53] 3个一维区间。为实现该范围查询，GS-Phoenix构造如下SQL语句：

图  3  6位Geohash编码

Figure 3.  6 Bits Geohash Encoding

SELECT c1, c2 FROM mytable WHERE geohash$\geq{}$16 AND geohash$\leq{}$31 UNION ALL

SELECT c1, c2 FROM mytable WHERE geohash$\geq{}$48 AND geohash$\leq{}$49 UNION ALL

SELECT c1, c2 FROM mytable WHERE geohash$\geq{}$52 AND geohash$\leq{}$53

该SQL语句包含3个子SQL语句，分别实现了3个Geohash一维区间子查询。GS-Phoenix将该SQL请求发送给Phoenix，最后将3个子查询结果合并到一起返回给客户端。

GS-Phoenix实现空间范围查询的过程如算法1所示。假定要查询矩形区域query_area内所有的空间点，首先计算左下角query_area.($x_{1}$, $y_{1}$)和右上角query_area.($x_{2}$, $y_{2}$)对应的两个Geohash编码，再计算这两个Geohash的最长公共前缀prefix。值得注意的是，prefix对应的空间矩形一定包含查询区域query_area。然后以query_area和prefix为输入参数，调用rect_to\_intervals算法（见算法2）获得查询区域query_area对应的多个Geohash取值区间并存放在链表intervals中。接着调用combine_intervals尽量将intervals中连续的区间合并成更大的区间。算法1针对intervals中的每个区间构造一个子SQL查询语句，最后将这些子SQL串拼接在一起，发送给Phoenix获得最终的查询结果。

算法2 rect_to_intervals是一个深度优先搜索算法。该算法的输入prefix对应一个矩形区域，且该区域可以用一个Geohash区间表示。在初始状态时，prefix对应的区域完全覆盖查询区域且面积大于查询区域，如果直接取prefix对应的Geohash区间作为过滤条件，往往会引入很多假阳性数据。因此不断地将输入参数prefix对应的矩形区域分解为对等的两个子区域，即分别为前缀prefix$\vert{}$$\vert{}$'0'和prefix$\vert{}$$\vert{}$'1'对应的区域。将该分解过程进行递归，直到子区域完全被包含在查询区域query_area内或者因子区域与query_area不相交被舍弃。此过程得到的每个子区域都对应一个Geohash区间。

如果rect_to_intervals对搜索深度不加以控制，则会产生大量琐碎的子区域，每个子区域内只包含很少的空间点或者不包含任何点。每个子区域对应Phoenix的一个子查询，最终将触发HBase的一个Scan操作。如果数据量太小，Scan操作的调度和启动时间将超过实际读取数据的时间，造成查询效率恶化。因此rect_to_intervals利用参数MAX_LENGTH限定不断增长的前缀串prefix的长度，这相当于限定了搜索深度，当被搜索的子区域小于一定面积时，将停止对其进一步分解。如果设定MAX_LENGTH为4，图 3中红色矩形所示的查询区域被分解为[16, 31]，[48, 51]和[52, 55] 3个一维区间，图 3中标深、浅阴影的两个子区域（[48, 51]和[52, 55]对应的区域）没有被进一步分解。

限定搜索深度可能会引入假阳性数据，如图 3区间[48, 51]对应的区域中50、51两个单元格的数据。为了消除假阳性数据，GS-Phoenix定义名为rect_contain(rect, lo, la)的UDF函数，该函数判断空间点(lo, la)是否位于rect定义的矩形区域内，并返回对应的布尔值。利用一维Geohash区间过滤条件粗选出部分空间点后，再利用UDF函数进一步过滤掉不符合条件的点。该精细过滤过程是在服务器端实现的，因此减轻了客户端的运算负载。对应可能会引入假阳性数据的Geohash区间，GS-Phoenix按如下方式构造子SQL查询语句：

SELECT c1, c2 FROM mytable WHERE geohash$\geq{}$48 AND geohash$\leq{}$51 AND rect_contain(query_area, lo, la)

为了进一步减小Scan操作调度和启动时间所占比例，在不引入假阳性数据的前提下，合并rect_to_intervals产生的部分Geohash区间，以减少Scan操作数量。如图 3虚线标识的查询区域在初始时被分解为[32, 35]、[36]、[38]、[40, 43]、[44]和[46]共6个区间，它们可以被合并为[32, 36]、[38]、[40, 44]和[46]共4个区间。算法1通过combine_intervals子程序实现了区间合并过程。

在进行空间范围查询时，用户还可以用任意多边形或圆形定义查询区域，也可以定义关于点、线、面等几何对象的缓冲区查询，即在点、线、面对象周围建立一定宽度的区域作为查询区域。GS-Phoenix处理不规则空间范围的基本方法是：首先求不规则查询区域的最小包围矩形(minimum bounding rectangle, MBR)，然后基于MBR执行§2.3的矩形范围查询进行粗过滤，最后利用与不规则查询区域类型匹配的UDF函数实施服务器端的细过滤。

GS-Phoenix实现空间kNN查询时主要利用了Phoenix的服务器端数据排序和数据聚集等优化机制。在处理包含ORDER BY和LIMIT n的SELECT查询时，Phoenix会在服务器端对粗选数据集进行排序并只保留前n条记录返回给客户端。如果该SELECT语句触发了多个Scan子查询或者涉及到多个Region服务器的子查询，Phoenix首先取每个子查询的前n条记录，然后在服务器端进行聚集，最后只返回聚集后的结果。假定要查询图 3中[48, 51]区间对应的区域内距离空间点Q最近的k个临近点，可以构造如下SQL语句：

SELECT c1, c2, distance(lo, la, Q) FROM mytable WHERE geohash$\geq{}$48 AND geohash$\leq{}$51 ORDER BY distance(lo, la, Q) LIMIT k

其中，distance(lo, la, Q)是用户自定义的UDF函数，用于计算空间点Q和点(lo, la)之间的欧氏距离。

考虑到Phoenix上述优化机制，GS-Phoenix实现kNN查询的基本方法是在以Q为中心点的合适的圆形区域内执行§2.4的空间范围查询，所不同的是构造SQL查询语句时包含"ORDER BY distance(lo, la, Q) LIMIT k"子句。该查询将返回已经按与Q点距离排序的至多k个空间点。如果查询得到的点数小于k，则扩大以Q为中心点的圆形区域大小，并作进一步的带排序的范围查询。GS-Phoenix实现kNN查询的过程如算法3所示。

算法3中，参数D等于估算的单位面积内的空间点数。如果数据集在空间上是均匀分布的，则D等于总点数除以空间总面积。算法3首先估算以Q为中心的圆形搜索区域的半径$r=\sqrt{(\Delta{}\cdot{}k)/(d\cdot{}\pi{})}$，使得在该圆内至少有k个点。参数$\Delta{}$是放大系数，用于适当扩大搜索面积，实验中取1.2或1.5。circle_to_intervals子程序用于将圆心为Q、半径为r的圆形搜索区域转化为多个一维Geohash区间。子程序knnquery_by_intervals则以这些Geohash区间为限定条件，以circle_contain UDF函数为细过滤条件，利用Phoenix的服务器端排序和聚集机制实现kNN查询。如果当前圆内的空间点数(等于len(result_points))小于k，则重新估算该区域内的点密度参数d，并重新计算圆形搜索区域的半径r，然后在新的圆内执行kNN查询。重新估算点密度d时，主要参考了当前圆形搜索区域内实际存在的点数，即len(result_points)的值，该点数更加真实地反映了Q点附近区域的点密度分布情况。

例如在图 3所示的空间内，中心点Q位于标号为12的单元格内。首先在以Q为中心的内层小圆内执行kNN查询，如果获得的空间点数小于k，则将搜索区域扩大到外层大圆。在对外层大圆进行kNN查询时，有些Geohash区间之前已经被搜索过，因此算法3通过difference子程序将重叠的区间删除后再进行搜索。

实际应用中，数据集在空间分布上往往是不均匀的。如果能够获得数据集在空间上的实际分布情况，那么就可以对算法的查询区域分解过程作进一步优化。假定子空间包含M条记录时，Scan查询调度启动的代价和假阳性数据细过滤的代价可以达到均衡。参数M类似于R树或QR树索引中控制节点分裂的参数。当算法2不断对搜索区域进行二分时，如果子区域内包含的点数小于等于M，则停止分解；否则，分解过程继续。因此，优化算法中只需要将算法2的第6步改成"IF estimate_num(prefix) $\leq{}$ M THEN"。其中，estimate_num(prefix)子程序的功能是根据数据集的空间分布统计估算prefix前缀对应区域内的空间点数。同理，在实现kNN查询时，基于实际空间分布统计便可以更加准确地估算圆形搜索区域的半径。

GS-Phoenix采用子表分裂的统计方法估算数据集的空间分布。HBase作为GS-Phoenix的底层数据存储系统，支持数据表自动分区机制，即当数据表中的数据超过一定长度时（默认为128 MB），则将该表拆分成两个子表。数据表中的记录是按照Geohash值由小到大排序的，根据Geohash编码的空间临近性可知，同一个二维子空间内的点在HBase表中也是临近存储的。如果某区域内的点密度大，则该区域对应的数据集会被自动拆分到多个子表。反之，如果某区域内的点密度小，则很可能将该区域及临近区域内的数据压缩到同一个子表中。通过HBase的客户端API，很容易获得一个HBase表的子表分裂情况以及每个子表的起始主键值（start_key）、终止主键值（end_key）和该子表内的记录数目n。通过计算Geohash区间[start_key, end_key]对应的子空间集合，假定n条记录均匀分布在这些子空间内，便可以获得这些子空间的密度统计数据。综合各个子表的统计数据，可以获得数据集所覆盖的整个空间的密度统计情况。

本文通过多项实验对GS-Phoenix的数据插入和空间查询效率进行评估。在搭建Hadoop及Phoenix分布式平台时采用4台Dell 7910 Tower工作站，每台配置双路Intel Xeon E5-2620 CPU（共32个虚拟核）和64 GB内存。每台工作站利用XenServer虚拟化软件创建多个虚拟计算节点。为了尽量减少磁盘竞争，让每个虚拟节点都工作在独立物理磁盘上。实验Hadoop平台中每个节点的主要软件配置为：CentOS6.8(64 bit)操作系统，Hadoop2.7.3和Phoenix4.7.0。采用仿真实验生成了两类带空间属性的数据，其中一类在空间上均匀分布；另一类在空间上符合Zipf分布^[14]，用于模拟具有多个热点区域的非均匀空间分布。

为了测试GS-Phoenix的数据插入效率，分别在具有4、8和16个数据节点的GS-Phoenix系统上测试插入n条数据记录所用的延时，并令n以10$^{7}$的速度递增。每种配置环境都测试了将Geohash列作为主键和作为二次索引的两种方案。客户端启动10个线程，同时向GS-Phoenix插入空间数据。每个线程负责插入n/10的数据，并记录各自完成任务所用的时间。各个线程记录的延时中，最长的延时作为插入n条记录的总延时。图 4给出了不同条件下插入n条数据所用的时间。

图  4  数据插入延时随记录数目n的变化情况

Figure 4.  Changes of Data Elapsed Time when Inserting n Records

由图 4可以看出，在每种条件下，数据插入延时随着记录数n的增大几乎都是呈线性增长的，说明GS-Phoenix系统具有比较稳定的插入吞吐率。随着数据节点数的增加，系统的插入吞吐率也会持续提高，这是因为插入数据的工作负载会被更多的数据节点共同分担。在16个节点的环境中，数据插入吞吐率达到最高，约17万/s，该吞吐率可以满足实时产生位置点数据的应用场景。

在进行空间范围查询实验时，把数据集中的点坐标都归一化到边长为1的正方形空间内。§2.3提出将查询区域分解为多个Geohash段进行粗过滤并配合UDF函数进行细过滤，为了验证该方法的查询效率，采用的对比方法是将查询区域映射为单个Geohash段进行过滤的直观方法。实验主要测试了0.01‰、0.1‰、1‰和5‰共4种不同选择率的空间范围查询的平均响应时间。为了准确对应选择率，将空间查询范围分别设成面积为0.003 2$^{2}$、0.01$^{2}$、0.032$^{2}$、0.071$^{2}$的正方形，并采用均衡分布的具有5 000万条记录的数据集。实验用的GS-Phoenix系统由16个节点构成。

图 5给出了空间范围查询平均响应时间的测试结果。其中，简单Geohash过滤和多重Geohash过滤的查询方案都采用了在客户端进行二次过滤的方法，而多重Geohash过滤+UDF过滤则利用UDF函数在服务器端实施二次过滤。通过实验结果可以看出，将空间查询区域映射为多个Geohash区间进行过滤的方法由于在粗过滤阶段去除了更多的冗余数据，所以其查询效率远高于简单Geohash过滤，尤其是在选择率较低(图 5中选择率为5‰的情况)且查询结果数据集较大时，其效率优势更加明显。图 5还验证了借助UDF函数在服务器端进行二次过滤可以进一步缩短查询响应时间，其原因是减少了假阳性数据的网络传播开销，且降低了客户端的计算压力。

图  5  空间范围查询平均响应时间

Figure 5.  Average Response Time for Spatial Range Queries

基于空间分布统计的查询优化方法实验中，在[0, 1]$^{2}$归一化空间中仿真生成了符合Zipf分布（参数s=1.1）的5 000万个空间点数据，并模拟了两个热点，其空间原始分布如图 6(a)所示。将测试数据均分成500个区，仅通过每个分区的起始行键和终止行键估算数据的空间分布，估算结果如图 6(b)所示。通过图 6可以看出，客户端简单的分区统计方法也能在很大程度上获得原始数据的空间分布情况。

图  6  原始空间分布与统计得到的空间分布

Figure 6.  Original Spatial Distribution and Statistical Spatial Distribution

针对该仿真数据集进行空间范围查询测试。对于0.01$^{2}$、0.02$^{2}$、0.03$^{2}$、0.05$^{2}$和0.07$^{2}$共5种大小的方形区域，分别采用未考虑空间分布和考虑空间分布的优化查询算法进行空间范围查询，对比平均查询响应时间，如图 7所示。由图 7可以看出，无论对于哪种大小的查询区域，根据所在区域的数据分布密度划分Geohash段的优化查询算法在效率上均优于不考虑空间分布的简单查询处理。

图  7  优化前后平均查询响应时间

Figure 7.  Average Response Time for Range Queries Before and After Optimizing

针对§3.3生成的非均匀分布数据集进行kNN查询实验。本实验测试了k分别等于1、10、100、1 000和10 000，参数$\Delta{}$分别等于1.2、1.5和2时，实施1 000次kNN查询的平均响应时间（图 8）。k=1时，等价于空间点查询。根据图 8的实验结果可知，对于常见的k值，其kNN查询都能在1 s内完成，因此GS-Phoenix可以满足实时查询和交互式数据分析的应用需求。

图  8  kNN查询平均响应时间

Figure 8.  Average Response Time for kNN Queries

本文设计了一种分布式海量地理空间数据管理系统GS-Phoenix，该系统构建在两种流行的开源项目Phoenix和HBase之上，其既具有NoSQL数据库的高可扩展、强容错、高吞吐率等特点，又支持带空间扩展的SQL语言对数据进行灵活访问。GS-Phoenix能够满足多类空间位置相关应用中关于实时插入数据和实时查询空间数据的需求。下一步将增加GS-Phoenix系统对空间线和面数据的存储和查询支持，另外还可以通过对Phoenix进行深度扩展，使其二次索引支持R-Tree或QR-Tree等更高效的空间索引方式。