CN113128607B - 一种轨迹相似度计算方法 - Google Patents

Abstract

本发明公开了一种轨迹相似度计算方法，本发明为了度量轨迹间的时空相似度，提出了一个新颖的RSTS模型来学习轨迹代表。在模型中考虑了轨迹中的时间成份。采用提出的时空感知的损失函数后，该模型可以学习到暗藏在大量历史轨迹中的精确移动路线中的位置转移模式，并且将这些学习到的轨迹时空特征编码到轨迹代表向量中用于轨迹相似度计算。通过开展了广泛的实验研究，实验结果证实了基于本方法学习到的轨迹代表进行轨迹相似度计算，对低质量的轨迹仍然具有鲁棒性。

Description

一种轨迹相似度计算方法

技术领域

本发明涉及数据挖掘技术领域，具体涉及一种轨迹相似度计算方法。

背景技术

随着位置追踪设备(如车辆导航系统和智能手机)和基于GPS的应用(如谷歌地图，百度地图，滴滴)的持续发展，轨迹数据的数量正在飞速增长。轨迹相似度度量作为时空数据库分析中的一个基本功能，已经被大量学者们广泛研究。面向多种不同或相同的应用场景，大量的方法被提出用于度量轨迹间的相似度。

已经存在的方法通常假设采样获得的轨迹具有一个统一的、连续的采样率，如果两条轨迹中的大部分采样的位置点都能形成一对一的点匹配，那么认为这两条轨迹是相似的。然而，对于不同的采样设备，由于各种各样的原因，采样率往往会发生变化。导致采样率不统一的原因，包括但不限于电池容量限制，间歇性的信号中断，以及不同的系统设定。在这种情况下，基于点匹配的方法是无效的，EDwP方法因此被提出通过动态插入位置点来进行轨迹匹配，以此来解决位置点采样率不统一的问题。接着，t2vec借鉴深度学习的方法，通过学习轨迹代表来度量轨迹间的相似度。然而，在t2vec模型中，轨迹中的时间信息被忽略了，由于没有衡量轨迹间时间维度的差异，因此它不能用于检索时空数据库中基于时空相似度的查询。在这种查询中，轨迹中的时间信息和空间信息显得同样重要，当计算轨迹的时空相似度时，这两种信息均是不可或缺的。

发明内容

针对现有技术中的上述不足，本发明提供的一种轨迹相似度计算方法解决了现有轨迹相似度计算方法未同时考虑时间信息和空间信息的问题。

为了达到上述发明目的，本发明采用的技术方案为：

提供一种轨迹相似度计算方法，其包括以下步骤：

S1、对原始轨迹数据进行预处理后通过添加噪声生成对应的噪声轨迹数据；建立初始的基于代表的时空相似度计算模型；基于代表的时空相似度计算模型包括编码器和解码器；编码器和解码器均包括若干依次相连的神经网络；

S2、将原始轨迹数据和其对应的噪声轨迹数据分布转化为时空格子编码序列，对应得到时空格子原始编码序列和时空格子噪声编码序列；

S3、将一条时空格子噪声编码序列作为当前的基于代表的时空相似度计算模型的编码器的输入进行前向计算，得到代表向量；

S4、将代表向量和时空格子原始编码序列作为当前的基于代表的时空相似度计算模型的解码器的输入进行前向计算，得到预测的序列和隐藏层状态序列；

S5、基于预测的序列、隐藏层状态序列和损失函数，采用梯度下降法更新当前的基于代表的时空相似度计算模型的参数；

S6、计算当前的基于代表的时空相似度计算模型的损失；

S7、判断损失是否降低，若是则保留当前的基于代表的时空相似度计算模型并返回步骤S2；否则进入步骤S8；

S8、判断损失未降低的次数是否大于阈值，若是则输出当前的基于代表的时空相似度计算模型并进入步骤S9；否则返回步骤S2；

S9、采用当前的基于代表的时空相似度计算模型进行轨迹相似度计算。

进一步地，步骤S1中对轨迹数据进行预处理的具体方法为：

将长度不在设定范围的原始轨迹数据进行删除，将空间坐标不在设定范围的原始轨迹数据进行删除，完成对原始轨迹数据的预处理。

进一步地，步骤S1中通过添加噪声生成对应的噪声轨迹数据的具体方法为：

对于预处理后的原始轨迹数据，采用5种不同的下采样概率生成5条下采样轨迹；对于每一条下采样轨迹，采用4种不同的噪声偏移概率添加噪声生成4条噪声轨迹，即每一条原始轨迹数据对应20条噪声轨迹。

进一步地，步骤S2中将轨迹数据转化为时空格子编码序列的具体方法包括以下子步骤：

S2-1、以地图的经度和维度作为X轴和Y轴、以时间作为Z轴构建三维空间，将三维空间划分为若干同样大小的格子，得到时空格子，并对时空格子进行编号；

S2-2、根据轨迹数据中的位置点和时间信息将该轨迹数据中的每个轨迹点与时空格子进行对应，将时空格子的编号作为与其对应的轨迹数据的编码值；

S2-3、将出现频次小于阈值的轨迹数据编码值替换为其邻居编码值，对所有编码值重新计数并按照从小到大重新编号；

S2-4、将选取的轨迹数据补齐至该批轨迹数据中最长轨迹的长度，并在当前每条轨迹数据的首位添加一个首位标志；

S2-5、将当前每条轨迹数据的位置点替换成与其对应的最新编号值，并在该轨迹数据的末尾添加末尾标志，得到时空格子编码序列；其中所有最新编号值构成词汇表。

进一步地，步骤S2-4中补齐的具体方法为添加负数；步骤S2-4中的首位标志和步骤S2-5中的末尾标志均为负数；其中补齐用的负数与首位标志用的负数不同，补齐用的负数与末尾标志用的负数不同。

进一步地，步骤S3中基于代表的时空相似度计算模型的编码器进行前向计算的具体方法为：

将输入的编码序列的第一个编码值输入编码器的第一个神经网络B₁，将输入的编码序列的第n个编码值和第n-1个神经网络B_n-1的隐藏层状态输入编码器的第n个神经网络B_n中获取编码器第n个神经网络B_n的隐藏层状态，将编码器的最后一个神经网络的输出作为代表向量；其中n≥2。

进一步地，步骤S4中基于代表的时空相似度计算模型的解码器进行前向计算的具体方法为：

将代表向量和时空格子原始编码序列中的首位标志输入解码器的第一个神经网络J₁中，将解码器的第m-1个神经网络J_m-1的隐藏层状态和时空格子原始编码序列中第m-1个编码值输入解码器的第m个神经网络J_m中获取解码器的第m个神经网络J_m的隐藏层状态h_m和预测值进而得到预测的序列/>和隐藏层状态序列H＝{h₁,h₂,...,h_m,...,h_M}；其中m大于等于2，M表示解码器的神经网络总数，/>和h_M分别为解码器的第M个神经网络J_M的预测值和隐藏层状态。

进一步地，步骤S5中损失函数L为；

L＝L₁+β×(L₂(a₁,p₁,n₁)+L₂(a₂,p₂,n₂))

dist(u,y_t)＝(1-λ)dist_s(u,y_t)+λdist_t(u,y_t)

L₂(a₁,p₁,n₁)＝max{d(a₁,p₁)-d(a₁,n₁)+margin,0}

L₂(a₂,p₂,n₂)＝max{d(a₂,p₂)-d(a₂,n₂)+margin,0}

其中y_t为需要预测的目标轨迹点，N_st为与y_t在时间和空间上最接近的邻居集合；u为编码值，即常数；h_t为解码器的第t个神经网络的隐藏层状态；V为词汇表；W^T为将隐藏层状态从隐藏层空间映射到词汇表空间的映射矩阵；表示W^T中的第u行；dist(·)表示dist函数，dist_s(·)表示空间距离，dist_t(·)表示时间距离；θ为权重值，θ∈(0,1)；λ为权重值，λ∈[0,1]；y为需要预测的目标轨迹的长度；u∈N_st表示编码值u为y_t在时间和空间上最接近的邻居中的个体；d(·)表示距离；a₁和p₁为与同一条时空格子原始编码序列对应的两条不同的时空格子噪声编码序列，n₁为其他任一条时空格子原始编码序列对应的一条时空格子噪声编码序列，且n₁与a₁在时空上的距离大于p₁与a₁在时空上的距离；a₂、p₂和n₂均为不同时空格子原始编码序列对应的时空格子噪声编码序列，且n₂与a₂在时空上的距离大于p₂与a₂在时空上的距离；margin为常数；β为权重值，β∈(0,1)。

进一步地，步骤S6的具体方法包括以下子步骤：

S6-1、构建测试集；测试集中包括Q个时空格子原始编码序列，以及与Q个时空格子原始编码序列对应的时空格子噪声编码序列；

S6-2、分别将测试集中的时空格子原始编码序列和时空格子噪声编码序列作为当前的基于代表的时空相似度计算模型的输入，对应得到若干组代表向量v₁和代表向量v₂；

S6-3、根据公式：

计算代表向量v₁和代表向量v₂的KL散度损失，将测试集对应的所有KL散度损失的平均值作为当前的基于代表的时空相似度计算模型的损失；其中I表示代表向量的维度；v₁(x_i)表示代表向量v₁中的第i个元素；v₂(x_i)表示代表向量v₂中的第i个元素。

进一步地，步骤S9的具体方法包括以下子步骤：

S9-1、采用与步骤S2相同的方法分别为给定的两条轨迹编码一个时空格子编码序列；

S9-2、将给定的两条轨迹对应的时空格子编码序列输入当前的基于代表的时空相似度计算模型并得到对应的两条代表向量；

S9-3、将给定的两条轨迹对应的两条代表向量的欧式距离的倒数作为轨迹相似度，即欧氏距离越大，轨迹相似度越低，完成轨迹相似度计算。

本发明的有益效果为：本方法通过把空间和时间分别划分为相同大小的空间格子和相同大小的时间片，先在二维空间上对空间格子进行编号，再在时间片上进一步编号，从而每一条轨迹就转化为了一个离散的数值序列；接着利用一个时空感知的损失函数结合三元损失来共同训练模型，当要解码一个目标位置点时，该时空感知的损失函数倾向于鼓励解码器分配一个更高的可能性给离该目标位置点时间和空间上最近的邻居位置点，模型以大量的历史轨迹作为输入，通过训练，精确移动路线中隐藏的位置转换模式便可以被有效学习到，进而在保留时间和空间信息的基础上实现轨迹相似度计算。

附图说明

图1为本发明的流程示意图；

图2为变化轨迹P的数目大小、下采样率和噪声偏移率时的最相似轨迹搜索结果示意图；

图3为当K＝100，200，300时，变化下采样率和噪声偏移率时k-nn准确率示意图；

图4为三元损失函数效果评估示意图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，该轨迹相似度计算方法包括以下步骤：

S1、对原始轨迹数据进行预处理后通过添加噪声生成对应的噪声轨迹数据；建立初始的基于代表的时空相似度计算模型；基于代表的时空相似度计算模型(RSTS模型)包括编码器和解码器；编码器和解码器均包括若干依次相连的神经网络；

S2、将原始轨迹数据和其对应的噪声轨迹数据分布转化为时空格子编码序列，对应得到时空格子原始编码序列和时空格子噪声编码序列；

S3、将一条时空格子噪声编码序列作为当前的基于代表的时空相似度计算模型的编码器的输入进行前向计算，得到代表向量；

S4、将代表向量和时空格子原始编码序列作为当前的基于代表的时空相似度计算模型的解码器的输入进行前向计算，得到预测的序列和隐藏层状态序列；

S5、基于预测的序列、隐藏层状态序列和损失函数，采用梯度下降法更新当前的基于代表的时空相似度计算模型的参数；

S6、计算当前的基于代表的时空相似度计算模型的损失；

S7、判断损失是否降低，若是则保留当前的基于代表的时空相似度计算模型并返回步骤S2；否则进入步骤S8；

S8、判断损失未降低的次数是否大于阈值，若是则输出当前的基于代表的时空相似度计算模型并进入步骤S9；否则返回步骤S2；

S9、采用当前的基于代表的时空相似度计算模型进行轨迹相似度计算。

步骤S1中对轨迹数据进行预处理的具体方法为：将长度不在设定范围的原始轨迹数据进行删除，将空间坐标不在设定范围的原始轨迹数据进行删除，完成对原始轨迹数据的预处理。

步骤S1中通过添加噪声生成对应的噪声轨迹数据的具体方法为：对于预处理后的原始轨迹数据，采用5种不同的下采样概率生成5条下采样轨迹；对于每一条下采样轨迹，采用4种不同的噪声偏移概率添加噪声生成4条噪声轨迹，即每一条原始轨迹数据对应20条噪声轨迹。

步骤S2中将轨迹数据转化为时空格子编码序列的具体方法包括以下子步骤：

S2-1、以地图的经度和维度作为X轴和Y轴、以时间作为Z轴构建三维空间，将三维空间划分为若干同样大小的格子，得到时空格子，并对时空格子进行编号；

S2-2、根据轨迹数据中的位置点和时间信息将该轨迹数据中的每个轨迹点与时空格子进行对应，将时空格子的编号作为与其对应的轨迹数据的编码值；

S2-3、将出现频次小于阈值的轨迹数据编码值替换为其邻居编码值，对所有编码值重新计数并按照从小到大重新编号；

S2-4、将选取的轨迹数据补齐至该批轨迹数据中最长轨迹的长度，并在当前每条轨迹数据的首位添加一个首位标志；

S2-5、将当前每条轨迹数据的位置点替换成与其对应的最新编号值，并在该轨迹数据的末尾添加末尾标志，得到时空格子编码序列；其中所有最新编号值构成词汇表。

步骤S2-4中补齐的具体方法为添加负数；步骤S2-4中的首位标志和步骤S2-5中的末尾标志均为负数；其中补齐用的负数与首位标志用的负数不同，补齐用的负数与末尾标志用的负数不同。

步骤S3中基于代表的时空相似度计算模型的编码器进行前向计算的具体方法为：将输入的编码序列的第一个编码值输入编码器的第一个神经网络B₁，将输入的编码序列的第n个编码值和第n-1个神经网络B_n-1的隐藏层状态输入编码器的第n个神经网络B_n中获取编码器第n个神经网络B_n的隐藏层状态，将编码器的最后一个神经网络的输出作为代表向量；其中n≥2。

步骤S4中基于代表的时空相似度计算模型的解码器进行前向计算的具体方法为：将代表向量和时空格子原始编码序列中的首位标志输入解码器的第一个神经网络J₁中，将解码器的第m-1个神经网络J_m-1的隐藏层状态和时空格子原始编码序列中第m-1个编码值输入解码器的第m个神经网络J_m中获取解码器的第m个神经网络J_m的隐藏层状态h_m和预测值进而得到预测的序列/>和隐藏层状态序列H＝{h₁,h₂,...,h_m,...,h_M}；其中m大于等于2，M表示解码器的神经网络总数，/>和h_M分别为解码器的第M个神经网络J_M的预测值和隐藏层状态。

步骤S5中损失函数L为；

L＝L₁+β×(L₂(a₁,p₁,n₁)+L₂(a₂,p₂,n₂))

dist(u,y_t)＝(1-λ)dist_s(u,y_t)+λdist_t(u,y_t)

L₂(a₁,p₁,n₁)＝max{d(a₁,p₁)-d(a₁,n₁)+margin,0}

L₂(a₂,p₂,n₂)＝max{d(a₂,p₂)-d(a₂,n₂)+margin,0}

其中y_t为需要预测的目标轨迹点，N_st为与y_t在时间和空间上最接近的邻居集合，即一批轨迹位置点；u为编码值，即常数；h_t为解码器的第t个神经网络的隐藏层状态；V为词汇表；W^T为将隐藏层状态从隐藏层空间映射到词汇表空间的映射矩阵；表示W^T中的第u行；dist(·)表示dist函数，dist_s(·)表示空间距离，dist_t(·)表示时间距离；θ为权重值，θ∈(0,1)；λ为权重值，λ∈[0,1]；|y|为需要预测的目标轨迹的长度；u∈N_st表示编码值u为y_t在时间和空间上最接近的邻居中的个体；d(·)表示距离；a₁和p₁为与同一条时空格子原始编码序列对应的两条不同的时空格子噪声编码序列，n₁为其他任一条时空格子原始编码序列对应的一条时空格子噪声编码序列，且n₁与a₁在时空上的距离大于p₁与a₁在时空上的距离；a₂、p₂和n₂均为不同时空格子原始编码序列对应的时空格子噪声编码序列，且n₂与a₂在时空上的距离大于p₂与a₂在时空上的距离；m arg in为常数；β为权重值，β∈(0,1)。

步骤S6的具体方法包括以下子步骤：

S6-1、构建测试集；测试集中包括Q个时空格子原始编码序列，以及与Q个时空格子原始编码序列对应的时空格子噪声编码序列；训练集和测试集可以按照8：2的比例进行分配；

S6-2、分别将测试集中的时空格子原始编码序列和时空格子噪声编码序列作为当前的基于代表的时空相似度计算模型的输入，对应得到若干组代表向量v₁和代表向量v₂；

S6-3、根据公式：

计算代表向量v₁和代表向量v₂的KL散度损失，将测试集对应的所有KL散度损失的平均值作为当前的基于代表的时空相似度计算模型的损失；其中I表示代表向量的维度；v₁(x_i)表示代表向量v₁中的第i个元素；v₂(x_i)表示代表向量v₂中的第i个元素。

步骤S9的具体方法包括以下子步骤：

S9-1、采用与步骤S2相同的方法分别为给定的两条轨迹编码一个时空格子编码序列；

S9-2、将给定的两条轨迹对应的时空格子编码序列输入当前的基于代表的时空相似度计算模型并得到对应的两条代表向量；

S9-3、将给定的两条轨迹对应的两条代表向量的欧式距离的倒数作为轨迹相似度，即欧氏距离越大，轨迹相似度越低，完成轨迹相似度计算。

在具体实施时，一个移动物体的精确移动路线R＝{s₁，s₂，...，s_n}记录了该物体经过的连续位置点的时空信息s_i＝([p₁，...，p_m]，t)。每个经过的位置点s_i包含一个m维的向量[p₁，...，p_m]代表该位置点的空间特征(比如，经纬度信息，三维坐标等)，和一个时间戳t代表经过该位置点的时间。实际上，由于位置追踪设备不会连续的记录移动物体经过的每一个位置，精确的移动路线往往是不可获得的。一条轨迹τ＝{s₁，...，s_n}是一条有限长度的时间序列，它是从一个移动物体的精确移动路线采样得到的。为了方便处理，用S_i＝([x，y]，t)来代表一个在二维平面上的离散位置点，其中[x，y]代表空间坐标，而t代表经过该位置点时对应的时间戳。

给定一批轨迹数据T＝{τ₁，...，τ_n}，我们旨在学习它们的轨迹代表V＝{v₁，...，v_n}用作轨迹间强健的时空相似度计算。具体来说，基于代表的轨迹相似度计算需要满足以下两条限制：

if sim_s(τ_i，τ_j)≥sim_s(τ_i，τ_k)and sim_t(τ_i，τ_j)≥sim_t(τ_i，τ_k)，then dist(v_i，v_j)≤dist(v_i，v_k)；

if dist(v_i，v_j)≤dist(v_i，v_k)，then dist(v_i′，v_j′)≤dist(v_i′，v_k′)。

其中，sim_s(·)和sim_i(·)分别代表两条轨迹间的空间相似度和时间相似度，dist(·)代表两个向量间的距离如欧式距离。v是原始的采样轨迹τ的轨迹代表，而v′是含有采样损失和噪声偏移的轨迹τ′的轨迹代表。

RSTS模型的编码器和解码器都不是固定的，可选的有CNN/RNN/BiRNN/GRU/LSTM等等，在RSTS模型中，可以选择使用GRU作为RNN单元。

在本发明的一个实施例中，使用了两个真实的轨迹数据集，第一个轨迹数据集提取自北京出租车数据，它记录了北京城区内大约10000辆出租车一周的活动轨迹。第二个轨迹集来自葡萄牙波尔图，记录了超过19个月总共170万条车辆轨迹数据。在北京出租车数据集中，平均的轨迹样本采样间隔为177秒，采样点之间的平均距离为623米。在波尔图数据集中，每辆出租车每隔15秒便会回传自己的位置数据。为了在两个数据集中构造低质量的轨迹作为模型输入，我们使用下采样和添加噪声偏移的方法，以不同的下采样概率和噪声偏移概率，随机的删除轨迹中的某些位置点或者在时间维度和空间坐标上添加一定的高斯噪声。为了保证下采样和添加噪声偏移后的轨迹仍然是有效的，我们规定原始轨迹的长度需要在20-100之间，同时经纬度坐标也需要在一定的范围之内。此外，我们将记录的时间戳限制在24小时之内。对于这两个轨迹数据集，我们都将前80％的轨迹用作训练集而将剩下的20％用作测试集。对于原始轨迹长度和经纬度的筛选条件见表1。

表1

	北京	波尔图
			经度范围	[116.250，116.550]	[-8.735，-8.156]
纬度范围	[39.830，40.030]	[40.953，41.307]
			长度范围	[20，100]	[20，100]

我们排除北京轨迹集中少于30个采样点经过的时空格子和波尔图轨迹集中低于50个采样点经过的时空格子，最终我们分别得到了15496和65718个经过次数较多的时空格子。每个采样的位置点是由离该采样点最近的这些经过次数较多的时空格子表示的。为了生成(τ_a，τ_b)对作为训练数据，我们对每一条高质量的采样轨迹τ_b进行了如下的下采样操作和添加噪声偏移操作。首先，我们以不同的下采样概率r₁＝[0，0.2，0.4，0.5，0.6]随机去除轨迹中的某些点来生成轨迹τ_b的5条子轨迹。接着我们以不同的噪声偏移概率r₂＝[0，0.2，0.4，0.6]对这些下采样轨迹分别进行经纬度坐标和时间维度上添加噪声处理。最终，对于每一条原始轨迹τ_b，我们得到了20组用于训练的轨迹组合(τ_a，τ_b)。

表2给出了默认的训练参数设定：使用Pytorch平台来训练我们的模型，装载一个Nvidia 1660GPU。如果验证集损失超过20000次都还没有减小，终止训练。所有的基准方法均使用Python编写，在win10平台上运行，使用一个AMD Ryzen 5CPU(3.6GHz)和32GB的内存。

表2

我们分别在测试集中随机选择了10000条轨迹记作Q，随机选择m条轨迹记作P。接下来我们会评估m取值对实验结果的影响。对于每一条轨迹τ_i∈Q，我们交替的取出该条轨迹中的采样点生成两条子轨迹D_Q和D_Q′。我们对P执行同样的操作得到子轨迹集合D_P和D_P′。接下来，对于每条轨迹τ_a∈D_Q我们检索它在D_Q′∪D_P′中的top-k个最相近的轨迹，并计算τ_a′的排名。直觉上，τ_a′应该排名很靠前因为τ_a′和τ_a都是由同一个轨迹交替采样生成的。

如图2(a)和图2(d)所示(横坐标为大小，纵坐标为平均秩)，当我们变化m时，对于EDR_t方法和EdwP_i方法，平均排名指标都有非常明显的增长趋势。而对于RSTS模型这种增长趋势明显减弱了。这个实验证明了RSTS能够更好地处理大规模的轨迹数据集。如图2(b)和图2(e)(横坐标为下采样率，纵坐标为平均秩)，固定|D_Q′∪D_P′|＝40K，当下采样率增大时所有方法的表现效果都下降了，但和其他两个方法相比，RSTS模型表现最好。图2(c)和图2(f)(横坐标为噪声偏移率，纵坐标为平均秩)展示了当我们变化噪声偏移率时对平均排名指标的影响。同样的，我们固定|D_Q′∪D_P′|＝40K。和EDR_t方法以及EDwP_t方法相比，很明显RSTS模型对噪声敏感度更低。即使当r2＝0.5，在北京数据集和波尔图数据集上，RSTS模型的表现都没有明显降低。

采用cross similarity偏差作为我们的评价指标，记作csd，计算如下：

其中τ_b和τ_b′是两条不同的原始轨迹，d(τ_b，τ_b′)可以看作是轨迹τ_b和τ_b′的路面实况，因此csd越小表明度量的距离越接近于地面实况。τ_a(r)和τ_a′(r)分别是轨迹τ_b和τ_b′的两个变形，它们是通过一定的下采样概率或者噪声偏移概率r产生的。我们从测试集中随机选择10000条轨迹对(τ_b，τ_b′)来计算它们的平均csd。由于空间限制，我们只展示在波尔图数据集上的实验结果。表3和表4分别展示了当我们变化下采样率r₁和噪声偏移率r₂下的平均csd。与EDR_t和EDwP_t相比，我们观察到RSTS模型总是取得更小的平均csd，这说明我们度量的时空相似度更符合路面实况。值得注意的是EDwP_t有时会产生一个更小的csd，这有可能是因为EDwP_t本身也考虑了如何更好地处理采样率不统一的轨迹。

表3

表4

r2	0.1	0.2	0.3	0.4	0.5
						EDRt	0.020	0.033	0.085	0.083	0.096
EdwPt	0.019	0.027	0.054	0.065	0.078
						RSTS	0.021	0.038	0.037	0.045	0.067

从测试集中分别选择了两组不同的轨迹数据集，一组包含1000个轨迹的轨迹集Q和一组包含了10000个轨迹的轨迹集DB。查询轨迹库Q和轨迹数据库DB可以看作是两组高质量的轨迹。对于每一个轨迹τ_i∈Q，我们查找数据库DB中离他最近的k个轨迹作为它的地面实况。接下来，我们对Q和DB以一个随机下采样概率r₁或者一个噪声偏移概率r₂来生成一对低质量的轨迹集Q′和DB′。我们对Q′和DB′进行上述同样的k-nn查询操作。最终，我们计算两组查询中共同的k-nn邻居比例作为k-nn准确率。从图3我们观察到，当我们增大r₁和r₂时，所有方法的准确率都下降了。当r₂从0.3变化到0.4时，EDR_t和EDwP_t准确率急速下降。此时很明显RSTS模型表现最佳。

用RSTS^#来代表没有使用三元损失函数训练的RSTS模型。如图4所示，当我们变化轨迹集P的大小时，RSTS^#模型产生的平均排名结果有明显的增长趋势，而RSTS模型产生的平均排名结果增势很小。因此，本方法提出的三元函数损失可以有效地提升平均排名这一指标，从而本方法的模型可以用于处理大规模的k-nn查询。

综上所述，本发明为了度量轨迹间的时空相似度，提出了一个新颖的RSTS模型来学习轨迹代表。在模型中考虑了轨迹中的时间成份。采用提出的时空感知的损失函数后，该模型可以学习到暗藏在大量历史轨迹中的精确移动路线中的位置转移模式，并且将这些学习到的轨迹时空特征编码到轨迹代表向量中用于轨迹相似度计算。通过开展了广泛的实验研究，实验结果证实了基于本方法学习到的轨迹代表进行轨迹相似度计算，对低质量的轨迹仍然具有鲁棒性。

Claims (9)

1.一种轨迹相似度计算方法，其特征在于，包括以下步骤：

S1、对原始轨迹数据进行预处理后通过添加噪声生成对应的噪声轨迹数据；建立初始的基于代表的时空相似度计算模型；基于代表的时空相似度计算模型包括编码器和解码器；编码器和解码器均包括若干依次相连的神经网络；

S2、将原始轨迹数据和其对应的噪声轨迹数据分布转化为时空格子编码序列，对应得到时空格子原始编码序列和时空格子噪声编码序列；

S3、将一条时空格子噪声编码序列作为当前的基于代表的时空相似度计算模型的编码器的输入进行前向计算，得到代表向量；

S4、将代表向量和时空格子原始编码序列作为当前的基于代表的时空相似度计算模型的解码器的输入进行前向计算，得到预测的序列和隐藏层状态序列；

S5、基于预测的序列、隐藏层状态序列和损失函数，采用梯度下降法更新当前的基于代表的时空相似度计算模型的参数；

S6、计算当前的基于代表的时空相似度计算模型的损失；

S7、判断损失是否降低，若是则保留当前的基于代表的时空相似度计算模型并返回步骤S2；否则进入步骤S8；

S8、判断损失未降低的次数是否大于阈值，若是则输出当前的基于代表的时空相似度计算模型并进入步骤S9；否则返回步骤S2；

S9、采用当前的基于代表的时空相似度计算模型进行轨迹相似度计算；

步骤S2中将轨迹数据转化为时空格子编码序列的具体方法包括以下子步骤：

S2-1、以地图的经度和维度作为X轴和Y轴、以时间作为Z轴构建三维空间，将三维空间划分为若干同样大小的格子，得到时空格子，并对时空格子进行编号；

S2-2、根据轨迹数据中的位置点和时间信息将该轨迹数据中的每个轨迹点与时空格子进行对应，将时空格子的编号作为与其对应的轨迹数据的编码值；

S2-3、将出现频次小于阈值的轨迹数据编码值替换为其邻居编码值，对所有编码值重新计数并按照从小到大重新编号；

S2-4、将选取的轨迹数据补齐至该批轨迹数据中最长轨迹的长度，并在当前每条轨迹数据的首位添加一个首位标志；

S2-5、将当前每条轨迹数据的位置点替换成与其对应的最新编号值，并在该轨迹数据的末尾添加末尾标志，得到时空格子编码序列；其中所有最新编号值构成词汇表。

2.根据权利要求1所述的轨迹相似度计算方法，其特征在于，步骤S1中对轨迹数据进行预处理的具体方法为：

将长度不在设定范围的原始轨迹数据进行删除，将空间坐标不在设定范围的原始轨迹数据进行删除，完成对原始轨迹数据的预处理。

3.根据权利要求1所述的轨迹相似度计算方法，其特征在于，步骤S1中通过添加噪声生成对应的噪声轨迹数据的具体方法为：

对于预处理后的原始轨迹数据，采用5种不同的下采样概率生成5条下采样轨迹；对于每一条下采样轨迹，采用4种不同的噪声偏移概率添加噪声生成4条噪声轨迹，即每一条原始轨迹数据对应20条噪声轨迹。

4.根据权利要求1所述的轨迹相似度计算方法，其特征在于，步骤S2-4中补齐的具体方法为添加负数；步骤S2-4中的首位标志和步骤S2-5中的末尾标志均为负数；其中补齐用的负数与首位标志用的负数不同，补齐用的负数与末尾标志用的负数不同。

5.根据权利要求1所述的轨迹相似度计算方法，其特征在于，步骤S3中基于代表的时空相似度计算模型的编码器进行前向计算的具体方法为：

将输入的编码序列的第一个编码值输入编码器的第一个神经网络，将输入的编码序列的第n个编码值和第n-1个神经网络/>的隐藏层状态输入编码器的第n个神经网络/>中获取编码器第n个神经网络/>的隐藏层状态，将编码器的最后一个神经网络的输出作为代表向量；其中/>。

6.根据权利要求1所述的轨迹相似度计算方法，其特征在于，步骤S4中基于代表的时空相似度计算模型的解码器进行前向计算的具体方法为：

将代表向量和时空格子原始编码序列中的首位标志输入解码器的第一个神经网络中，将解码器的第m-1个神经网络/>的隐藏层状态和时空格子原始编码序列中第m-1个编码值输入解码器的第m个神经网络/>中获取解码器的第m个神经网络/>的隐藏层状态/>和预测值/>，进而得到预测的序列/>和隐藏层状态序列；其中m大于等于2，M表示解码器的神经网络总数，/>和/>分别为解码器的第M个神经网络/>的预测值和隐藏层状态。

7.根据权利要求1所述的轨迹相似度计算方法，其特征在于，步骤S5中损失函数L为；

其中为需要预测的目标轨迹点，/>为与/>在时间和空间上最接近的邻居集合；/>为编码值，即常数；/>为解码器的第t个神经网络的隐藏层状态；V为词汇表；/>为将隐藏层状态从隐藏层空间映射到词汇表空间的映射矩阵；/>表示/>中的第u行；/>表示dist函数，/>表示空间距离，/>表示时间距离；/>为权重值，/>；/>为权重值，；/>为需要预测的目标轨迹的长度；/>表示编码值u为/>在时间和空间上最接近的邻居中的个体；/>表示距离；/>和/>为与同一条时空格子原始编码序列对应的两条不同的时空格子噪声编码序列，/>为其他任一条时空格子原始编码序列对应的一条时空格子噪声编码序列，且/>与/>在时空上的距离大于/>与/>在时空上的距离；/>、/>和/>均为不同时空格子原始编码序列对应的时空格子噪声编码序列，且/>与/>在时空上的距离大于与/>在时空上的距离；/>为常数；/>为权重值，/>。

8.根据权利要求1所述的轨迹相似度计算方法，其特征在于，步骤S6的具体方法包括以下子步骤：

S6-1、构建测试集；测试集中包括Q个时空格子原始编码序列，以及与Q个时空格子原始编码序列对应的时空格子噪声编码序列；

S6-2、分别将测试集中的时空格子原始编码序列和时空格子噪声编码序列作为当前的基于代表的时空相似度计算模型的输入，对应得到若干组代表向量和代表向量/>；

S6-3、根据公式：

计算代表向量和代表向量/>的KL散度损失，将测试集对应的所有KL散度损失的平均值作为当前的基于代表的时空相似度计算模型的损失；其中I表示代表向量的维度；/>表示代表向量/>中的第i个元素；/>表示代表向量/>中的第i个元素。

9.根据权利要求1所述的轨迹相似度计算方法，其特征在于，步骤S9的具体方法包括以下子步骤：

S9-1、采用与步骤S2相同的方法分别为给定的两条轨迹编码一个时空格子编码序列；

S9-2、将给定的两条轨迹对应的时空格子编码序列输入当前的基于代表的时空相似度计算模型并得到对应的两条代表向量；

S9-3、将给定的两条轨迹对应的两条代表向量的欧式距离的倒数作为轨迹相似度，即欧氏距离越大，轨迹相似度越低，完成轨迹相似度计算。