0 引言

近来随着时空数据的价值在各个领域逐渐凸显，尤其在商业航空领域行业，客机使用全球定位系统(Global Positioning System, GPS)进行高频率的信息采集产生大量的数据；同时时空数据本身具有高维度性、轨迹连续性，使时空数据在实际应用面临查询效率上的挑战。

目前国内外在这个方面都有很多研究，基于B-tree，尤其是在R-tree^[1]和R^*-tree^[2]上的变形和发展最多。这是因为对较大规模数据进行索引，如采用非平衡的树形索引结构很难满足查询效率的要求，即使如此面对现有的大规模数据其依然在查询效率上存在很大的问题。TPR-Tree^[3]、TPR^*-Tree ^[4]和HTPR^*-tree^[5]就是基于两者的重要变形，它们致力于减小索引规模，同时提高索引的查找和维护效率，如TPR-Tree引入类似最小边界矩形(Minimum Bounding Rectangle, MBR) 的速度边界矩形(Velocity Bounding Rectangle, VBR) 的概念，用于减小索引的重叠区域进而提高查找效率，但VBR不断地移动会增加算法的复杂性，难以满足现在系统对实时性的要求。

同时面对海量时空数据，树形结构的索引方法存在大量闲置的内部节点，占用了大部分的索引空间，而且存在难以用并行化方式进行索引管理的问题^[6]。近些年，采用空间划分并编码以实现索引的思想逐渐受到重视并用于图像空间和地理空间的检索^[7-9]，其中文献[6]提出的D-Tree实质上是一个虚拟的无内部节点的树，它以一定的编码方式对包含数据的叶子节点进行编码。D-Tree避免了对内部节点的存储，同时也剪除了空的叶子节点，降低了存储成本，一定程度上提高了查询效率。然而D-Tree完全根据数据分布进行区域划分，面对频繁大量的数据插入时，每次都需在所有数据上重新建立索引，这在高频写入的实时系统中将耗费大量的时间和空间资源。

针对上述问题，本文提出了HDL-index (Hybrid spatio-temporal data indexing based on data and query-log distribution, HDL-index)索引方法。HDL-index对实时数据采用空间划分的思想，使数据平均分布在所划分的单元格内；同时有效利用历史查询记录，将在查询热点区域上进行查询的历史记录聚类后抽象出查询模型，然后根据模型的坐标位置和其查询粒度对整体查询区域进行分割从而建立索引。最终基于GeoHash编码进行两部分索引的合并^[10]，同时对数据采用在时间维度上的多层存储方式进行存储。

1 HDL-index整体思想

本文提出的索引机制结合实时数据和历史查询记录两个部分进行索引的建立。对于实时数据部分，HDL-index算法根据数据分布进行区域划分，尽量均衡划分得到的每个子区域所包含数据的个数，而且在区域划分的同时对划分得到的子区域进行GeoHash 编码。对于历史查询记录部分，HDL-index算法采用贪心策略抽象出历史查询记录的代表模型^[11]，然后根据这些模型的分布，将于之对应区域内默认粒度的索引单元格迭代四分到小于等于模型的粒度，本文称之为四分基准化。同样对得到的索引单元格进行GeoHash编码，最终将两部分索引合并作为新的索引待查，具体流程如图 1所示。

对于图 1中历史查询记录是指在系统应用的空间范围内进行区域查询的记录，而对于每一个查询记录对应的查询区域本文都将其用一个能覆盖它的最小矩形表示，以便实际应用中方便操作。

在空间数据的类型上，为了应对更加复杂的应用要求，例如对于运送危险化学物品的车辆，按照交通规定需要距离人群、住宅区、加油站一定距离行驶停放，这里就要求对整个车辆进行处理。本文将时空数据类型从点对象扩展到非零面积对象，对非零面积对象采用上面处理查询区域的方式以矩形表示。这时考虑到数据的一致性，在处理空间点数据对象时将它当作一个矩形处理。例如：任意一个空间点P_X=(x, y)，则它的矩形形式为R_X=(P_LL, P_RU), P_LL=P_RU=P_X(P_LL为左下角坐标，P_RU为右上角点坐标，下文同)，但在实际存储时只存储P_X，这样不但解决了数据结构不统一的问题，而且没有占用额外的存储空间。

按图 1中的流程最终合并得到的索引I将为下一批数据到来之前的查询提供索引。在数据存储方面，海量时空数据的存储压力大多来自数据密集区域。根据HDL-index建立的思想，在数据密集区域内，图 1中属于I_d索引的每个索引单元格的数据条数约等于按数据分布进行空间划分时设置的每个单元格所能包含数据量的最大阈值N_max;同时在存储每个索引单元格包含的数据之前，需采用2个长整型记录，分别记录索引编码值和数据总长度，以及用两个时间戳存储这些数据中的最小和最大的时间值共32个字节的标识信息，数据的实际存储形式如图 2所示。

为了对图 2给出的数据存储方式进行形式化的定义，假设每条数据大小为R_size字节，每个磁盘块大小为B_size字节，每个磁盘页有k个磁盘块，则可定义用于存储图 1中属于I_d的单个索引单元格包含数据所需的磁盘页个数D_page为：

${{D}_{\text{page}}}\approx \frac{{{N}_{\text{size}}}{{R}_{\text{size}}}+32}{k{{B}_{\text{size}}}}$

(1)

图 1中I_d和I_l两部分索引都是对相关区域进行4等分建立起来的，但由于数据并非绝对均匀分布，因而在I_d和I_l合并得到I时也就不可能将I_d中需要细化的单元格内的数据绝对4等分，所以可定义存储图 1中属于I的单个索引单元格内的数据所需磁盘页的个数D为：

$D\approx {{{D}_{\rm{page}}}}/{{{4}^{p}},p\in \bf{N}}\;$

(2)

在按图 2的进行数据存储时，对不同应用可以根据R_size设置B_size的大小;同时图 2所示的存储结构为每个批次的数据分配若干个独立的磁盘页，将数据按时间顺序进行连续存储，这样可有效地解决时间维度上的查询问题。

2 HDL-index的结构与算法 2.1 面向实时数据的索引建立

对于当前时间窗口的实时数据，本文采用空间划分的方式，将整个二维的数据区域沿两个维度的中线分成4个子矩形区域，然后按这种方法进行递归划分，最终使每个子矩形内的数据量尽量相当。为了使每个子矩形内的数据量相近，D-tree在进行区域划分时，需判断区域内数据的个数，使其小于等于给定的阈值N_max，但对于涉及多个对象的大规模实时数据系统，同一个坐标点在窗口时间内的不同时刻很可能产生多个数据。假设这个数值为n_r，当n_r≥N_max时递归将不能完成，为此需要添加一个参数L_min，限制子区域最小边长。另外，由于本文将数据类型推广到二维对象，当出现非零面积数据大面积交叠时递归同样不能完成，所以需要在划分中判断划分的四个子区域两两对角区域内的数据量较其父区域是否减少，如果两个对角区域都未减少，则说明存在重叠区域不必再继续划分。具体划分过程如算法1。

算法1 实时数据的划分(AreaDiv)。

输入：初始矩形R_init, N_max, L_min, 数据集合D, 初始矩形编码C_init。

输出：C_Res:Tuple t=(Code, Entry-Set)。

1) L_side ← float:GetRecSide(R_init);

2) subRecs ← QuadRec(R_init);

3) Define array D_sub record data set in the subRecs

4) exit ← NumberLE(D, subRecs, N_max, D_sub);

5) IF L_side≤L_minor exit is true THEN

6)  Add the data that contained in R_init into Entry-Set E_current;

7)  Instance a new t=(C_init, E_current), add it into Res;

8)  RETURN C_Res;

9) END IF

10) C_sub←C_init+"00";

11) AreaDiv (subRecs[0], N_max, L_min , D_sub[0], C_sub);

12) C_sub←C_init+"01";

13) AreaDiv (subRecs[1], N_max, L_min , D_sub[1], C_sub);

14) C_sub←C_init+"10";

15) AreaDiv (subRecs[2], N_max, L_min , D_sub[2], C_sub);

16) C_sub←C_init+"11";

17) AreaDiv (subRecs[3], N_max, L_min , D_sub[3], C_sub);

在算法1中，调用的函数NumberLE主要是进行区域内数据个数是否小于等于N_max的判断，用以解决上面提到的数据位置重叠的问题。Entry-Set中每个实体都代表一个空间数据对象，具体形式为(TraID, ObjPointer)，其中TraID为空间对象一次移动过程的轨迹编号，ObjPointer指向数据地址。算法1中，在区域划分的同时进行GeoHash编码，GeoHash的具体编码方式如图 3所示。

采用这样的方式进行编码，可以通过所提供的经纬度值快速定位到其所对应的索引编号，继而根据其他查询条件获得最终结果, 同时这样的编码方式对任意给定位置的邻近查询操作也是非常高效的。当所查找的相邻区域和给定位置处于同一级别的编码时，可直接获取索引编码；当两者不在同一级别时，由于GeoHash是前缀编码，所以依然能顺利得到所要的索引编码。

假设要查询图 3(b)中，与索引编码I_C=1100所包含的某个空间对象相距L_Q (单元格的边长为L_U，假设此时L_Q=L_U)的符合某些条件的空间对象。以查找I_C正下方索引为例，只需取I_C奇数位置上的值，与L_Q /L_U所得整部分值的二进制形式作二进制的减法运算，然后用所得结果替换I_C对应位置的值即可得到要查找的索引编码1001，对照图 3(b)可以看到1001正是I_C正下方索引编码。按照这样的计算方法，结合表 1给出的不同查找方位的操作列表即可求得其正上方、正左方和正右方的索引编码。

表 1中操作列中的值因二进制无法直接表述，因此采用的是十进制表述。在实际操作中需统一转化为二进制再进行运算。对于其他四个对角线方向的邻近查找可通过正方向间接计算出来，其中N=(2^n/2-1) ≤(n-1) ，确保运算可正确进行。

2.2 面向历史查询记录的索引建立

只面向数据通过空间划分的方式建立索引，在每次数据更新时都需要对所有有效数据重新进行划分。在大数据时代，面对频繁和大量的数据更新，如果采用上述的办法，将耗费很长的时间用于索引的建立;同时，只面向数据建立的索引只能体现数据的分布，并不一定能提供较好的查询效果，即数据的分布和查询的分布并不是总是一致的。在网络拓扑等相关研究中，利用查询记录辅助索引建立方面国内外学者已经有不少研究^[12-14]，同样的本文将历史查询记录分布作为构成索引的一部分，用以快速建立索引并提高查询效率。

对于时空数据的查询从空间角度可分为区域查询和点查询。本文将区域查询所覆盖的范围用其形心代表，这样可以将查询统一表示为元组〈P_Q, R_Q〉，其中P_Q为区域查询的形心或者查询点，R_Q为查询半径，对于区域查询设置其值为能覆盖查询区域最小正方形边长的1/2，对于点查询R_Q=0。在面向历史查询记录建立索引时，首先需要对所有P_Q采用DBSCAN(Density-Based Spatial Clustering of Applications with Noise)算法^[15]进行聚类，使得的每个cluster内每两个相邻的P_Q的距离不大于r(r为根据具体数据情况给定的经验值)，这样做的目的是为了使接下来抽象出查询模型能尽可能大地覆盖其代表的各个查询区域的面积，即保证模型尽量真实。然后在对每个cluster内的元组进行模型抽象时，采用改进的DBSCAN算法获取最近邻点的算法依次获取Q_model，它是以某个确定P_Q为圆心P_model，以r为半径圆的外接正方形，并且能尽可能多地包含P_Q。再用统计的方式获取Q_model中P_Q对应的出现比例为p_rate的R_Q值为R_rate。最后构造一个以P_model为中心，以L_c=2(R_rate+r)为边长的正方形作为查询代表模型，这样在给定p_rate的情况下，使抽象出的模型可以覆盖其所代表的绝大部分历史查询区域。这个具体过程如算法2所示。

算法2 抽象查询代表模型(CreateModel)。

输入：cluster Tuple t_c=〈P_Q, R_Q〉集合S, r, p_rate。

输出：M_Res:Tuple t_m=(Center-Point, Side-Length)。

1) IF the number of t_c in S is equal 0 THEN

2)  RETURN M_Res;

3) END IF

4) IF the number of〈P_Q, R_Q〉in S percent is equal p_rate THEN

5)  Instance a new t_m=(P_Q, 2(R_Q+r)) and add it into M_Res;

6)  RETURN M_Res:

7) END IF

8) Create a set S_r to keep rest of t_c after GreedyModel(S_r, r);

9) S_r ← S;

10) Instance a new t_m ← GreedyModel(S_r, r);

11) Add t_m into M_Res;

12) CreateModel(S_r, r);

算法2中函数GreedyModel(S_r, r)主要实现贪心策略以获取能代表最多历史查询的代表模型。HDL-index将所有抽象得到的查询代表模型记录在Json文件中，这样以配置文件的形式进行下面索引建立的操作，有助于索引的并行化。

在根据代表模型进行空间划分上，首先需要根据系统设置将整个查询区域按照2.1节的方法划分成若干区域，作为默认查询粒度为L_d的索引单元格。然后将所有的查询模型对应到查询空间上，检查各个查询代表模型的边长L_d与L_c的大小关系。如果L_d>L_c，则采用2.1节的方法将相关默认粒度的索引单元格进行递归划分，直至L_d≤L_c。由于划分的思想相似，只是面向历史查询记录的索引建立不涉及数据，只需要判断L_d和L_c的关系，不需要其他限制条件，在此就不再进行累述。最后将得到的结果同样采用GeoHash的编码方式进行编码，以便进行I_d和I_l两部分索引的合并。

2.3 索引的合并和更新策略

2.1和2.2节完成了图 1中I_d和I_l两部分索引的建立，按照HDL-index的整体思想，需要将两部分索引合并成最终索引。由于采用同样的编码方式为完成合并提供了前提，在合并过程中有以下三种情况：

1) 对应区域内I_d和I_l的编码相同，说明它们代表完全相同的查询区域，由于只有I_d对应实际数据，只需要保留I_d；

2) 对应区域内I_d的编码是I_l的编码的前缀，说明I_d的查询粒度是I_l的m倍(m为正整数)，这时需要将I_d的编码对应的数据，对应到以I_d的编码为前缀的I_l中相应的索引区域内并保存I_l；

3) 对应区域内I_l的编码是I_d的编码的前缀说明I_l查询粒度较大，保存I_d即可。

在合并中第2) 种情况要对I_d进行四分划分，划分后得到的索引数量必然呈4倍量级地增长，从而产生更多的索引编码，但是I_l代表查询的热点，所以此时小粒度的索引将大幅度提高查询效率。由于HDL-index具有较简单数据结构，在此是用较小空间换取了较高的时间效率。

经过上面的合并可以得到最终的索引，对于索引更新的基本策略如下：

1) 在系统初始化时，由于没有历史查询记录，HDL-index只采用I_d作为最终索引，等待第一批查询记录到来；

2) 当查询记录到来时，按照HDL-index的思想进行索引的建立得到完整索引，然后当实时窗口数据再次到来时，停止I_l的建立，只将I_d与完整索引合并；

3) 当步骤2) 不断合并过程中，上面合并操作第1) 种情况和第3) 种情况合并掉的编码总数量，等于步骤2) 中I_l的编码的数量时，回到步骤2) 。

对于HDL-index基本更新策略中步骤3) 触发条件实质上是指由于数据的不同和数据整体分布的改变使得查询热点变化，因而目前的索引已不能很好地适应当前的查询需求，所以需要对索引进行全部重建。由于最终索引结构与D-tree类似可采用相同的操作进行轨迹数据的删除，具体算法可参考文献[6]，不同的由于HDL-index采用时间分层存储，所以可以对无时效性的数据进行批量删除。

3 HDL-index性能评估与实验

为了评估HDL-index的性能，本文采用国内某航空公司百万兆字节级别的真实数据进行各项性能的测试。近数千万条的数据分布在我国中东部地区和东南亚地区，这些数据是该航空公司采用GPS以1s为频率采集的，实验是针对在E114°~E116°和N21°~N23°区域内的数据进行的，大致为深圳、香港及其沿海周边空运繁忙地区为重点实验区域，同时经分析发现这些数据是在以边长为1经纬度的区域内聚集分布的。

为了更清楚体现HDL-index性能的优缺点，本文选取上文提到的D-tree作为比较对象，因为两者都采用空间编码的方式建立索引。本文的实验在一台2.20 GHz Intel Core Duo CPU, 2 GB RAM计算上机进行，同时采用同样数据集和相关参数设置。

3.1 索引空间性能评估

为了进行占用空间大小评估，定义了下面一系列参数以方便计算。

M：占用空间总大小；

S：总数据量

L：划分得到的索引单元格个数；

C：索引单元格的最大容量；

p：一个指针的大小；

d：一个双精度浮点数的大小；

i：一个整数的大小。

对于HDL-index和D-tree虽然采用不同的编码方式，但由于都是对空间数据进行研究，两者具有相同的索引数据形式，每个索引单元格由一个(Code, Entry-Set)元组(其中Code整数编码；Entry-Set为数据实体集合)表示，Entry-Set中每个实体一个(TraID, ObjPointer)元组(其中TraID为整数组合编码，由设备号和唯一时间标识组成；ObjPointer为指针，指向数据)表示。所以两者空间占用大小M可以统一表示为：

$M=L*p*(C+1)+S*(p+8*d)$

(3)

在同样的实验数据下在表达式(3) 中，两者只有L和C两个参数不同，本文进一步将参数C取相同值，则决定两者占用空间大小只有参数L。由于D-tree直接采用Hash表的方式进行索引存储，而HDL-index采用时间分层的索引存储方式。实验将系统置于较大压力下运行，即每个批次的实时数据在划分后，都能使每个索引单元格达到其最大容量C，此时在k次数据到达后，则有L_D-tree=kL_HDL-index，图 4展示了数据批次数(每批约100000条数据)和索引空间的关系。

由于HDL-index是由I_d和I_l两个部分合并得到的，而合并后的索引单元格会增加，但是由于数据量远大于查询记录提取得到查询代表的个数，合并过程中索引单元格个数只会少量增加，所以L_D-tree的值会小于kL_HDL-index。

3.2 时间性能评估

在索引的建立上D-tree只从数据考虑方面，通过空间划分使得每个索引单元格内所包含的数据个数尽量均衡在N_max左右。按照这样策略进行索引的建立时，由于D-tree采用的是递归方法，所以D-tree耗费的时间t与数据量的个数是正相关的。HDL-index由于采用了在时间维度上的多层存储结构，使得HDL-index只需要对当前窗口数据进行划分。对于稳定系统的当前窗口数据，可以认为不同批次数据量是大致相同的，所以在HDL-index第一部分索引建立时间耗费上基本是一定的;同时由于HDL-index第二部索引只需读取Json配置文件，然后与默认查询粒度比较建立索引，这过程几乎不耗费时间。所以与D-tree相比，HDL-index建立索引的时间呈大致水平状态。图 5给出了数据量与索引建立耗费时间的关系的平均实验结果。

由图 5可以看到，只有在数据量较少时HDL-index两部分索引建立并合并的时间会大于D-tree。当数据量持续增大时，两者变化趋于稳定，但由于某些批次新增数据分布情况存在较大的变化对应地导致时间上的一定变化。

3.3 查询性能评估

由于HDL-index和D-tree都是采用编码思想建立索引的，所以在设置相同的N_max值时，两者在一般区域查询的效率基本一样。只需要通过查询条件计算出对应索引编码，再根据其他条件进行筛选，具体可参见文献[6]。但在热点区域查询方面，由于 D-tree只面向数据进行索引的建立，而HDL-index将历史查询记录考虑进来，在热点区域HDL-index根据查询将该区域进行更细粒度的划分，虽然会产生更多索引编码，但是对于时空应用绝大多数的查询都发生在热点区域，因此热点区域细粒度的索引不但大幅度提高该区域的查询效率，而且也提高了系统的平均查询效率。

对于热点区域划分得到索引单元格的个数和粒度，取决于抽象出的查询代表模型的位置和查询粒度(即模型的边长)，而它们分别受算法2中r和 p_rate两个参数影响。由本章开始所给出的数据分布和聚集特点，实验在所提到的相对较小的数据密集区域生成10000个窗口查询，这些查询的查询半径在0.05375~0.51557, 同时经分析发现数据聚集区域间隔在0.2经纬度左右，所以r的取值要小于0.2以保证每个聚类在同一个数据聚集区域内完成，具体r和p_rate对热点区域抽象模型的两方面影响如图 6所示。

由图 6中横向子图间对比可以发现，p_rate只影响模型的查询边长，当r值相同时，模型的位置和每个模型覆盖的查询个数相同(r=0.08时，模型个数为77；r=0.04时，模型个数为183) ;同时, 由图 6中纵向子图间对比可以发现r影响模型所覆盖的查询个数，此时在相同p_rate位置上的半径也因为其所覆盖的查询记录集合不同而不同，而本文使用贪心的方法去选取中心点，这种方法只关注所覆盖的查询数量，所以在实际实验结果中对应模型的位置也有稍微的偏移。

根据上面的分析，只要r<0.2时就能得到合理的模型，但是实验证明r值与最终得到的模型个数成反比，所以要根据应用数据设置r的大小，本文选取r=0.08和p_rate=0.75以进行下面的热点区域的查询实验，根据数据分布设置L_min=1，同时由于数据覆盖率较大结合实际划分发现将N_max设置为数据总量的8%左右时，能保证数据较好地平均到各个索引单元格内，然后进行不同热点区域查询实验得到平均查询时间，具体实验结果如图 7所示。由图 7可以看出，HDL-index的查询性能明显优于D-tree。这是由于选择图 6(b)的设置进行历史查询的模型生成，获得模型的查询粒度在0.1~0.92，而本文设置L_min=1，最终使得HDL-index在相关热点区域的查询粒度比D-tree精确了2~16倍，同时对于热点区域查询数据量越大HDL-index的查询性能优势明显。

4 结语

本文在已有的索引思想上，对大规模时空数据的索引进行了研究。HDL-index的提出将目前在研究时空索引上所忽略的查询记录，这一现成且具有实际意义的数据考虑进来，解决了已有时空索引方式存在的空间浪费和热点区域查询效率低的问题。本文经过理论分析和实验证明了HDL-index的可行性和有效性，但在处理大规模移动时空数据的轨迹相似性查询时还有较大不足。对此下一步准备引入超图对索引编码进行组织，以便能更好地处理这类查询，同时超图的引入也能更加清晰直观地展示不同轨迹之间的关系。

致谢: 在这里由衷地感谢清华大学数据科学研究院——工业大数据研究中心陆薇常务副主任提供的良好交流学习环境和实验设备，也特别感谢中心张硕博士在选题和研究过程中精心指导及建议。