阅读量:0
文章目录
前言
笔者在之前的研究中，尝试对车辆轨迹数据进行空间聚类，以期望发现车辆在行驶过程中的停留信息。在笔者之前的文章中，笔者使用了DBSCAN算法来做这一件事。
然而，对于时序的车辆经纬度数据，DBSCAN有一个很大的问题——没有考虑数据中蕴含的时间信息！时间信息是时间序列数据与其他数据区别的重要特征。举个例子：在使用DBSCAN对车辆的经纬度进行聚类时，它仅仅是把那些空间相近的数据聚成一类，可是同一类中的数据可能时间相差很大。我们想要发现车辆的驻留行为，那些空间和时间都相近的轨迹点才能聚成一类，这才表示这辆车可能在某段时间因为某原因发生了停留。
正因如此，很多年前的研究者就对DBSCAN进行改进，有了适合用作时间序列数据密度聚类的算法——ST-DBSCAN。STDBSCAN的具体算法本文不再赘述，基本流程和DBSCAN无异，只是在可达点寻找中加入了时间阈值作为限制条件。若想要了解详细算法，可自行互联网搜索。
本文还是以如下格式车辆轨迹数据为例，实提供了ST-DBSCAN对车辆轨迹数据聚类并分析的方法：
| collect_time | id | lon | lat |
|---|---|---|---|
| 时间 | 车辆标识 | 经度 | 纬度 |
为了尽量去除噪声影响，车辆轨迹数据已经经过滤波平滑，平滑方法可见作者之前文章：https://blog.csdn.net/jgsecurity/article/details/140608431。
一、单辆车轨迹的聚类与分析
为了尽量与scikit-learn库中的使用方法相似，本文用类来实现STDBSCAN。class STDBSCAN的内容可以放在单独文件中作为模块导入，也可以同一文件中使用。
1.引入库
使用了数学计算库numpy，数据分析库pandas，机器学习库scikit-learn，地理相关库shapely和geopy，绘图库matplotlib。
import numpy as np import pandas as pd from datetime import timedelta from shapely.geometry import MultiPoint from geopy.distance import great_circle from sklearn import metrics import matplotlib matplotlib.use('TkAgg') import matplotlib.pyplot as plt import matplotlib.colors as mcolors 2.class STDBSCAN实现
采用class来实现STDBSCAN。
类有四个属性：spatial_threshold（距离阈值，单位米）、temporal_threshold（时间阈值，单位分钟）、min_neighbors（邻域内最少点数）、labels_（聚类后的标签）。在__init__构造函数中为前三个参数设置了默认值。
retrieve_neighbors(self, index_center, df)用于寻找给定一个核心点的所有可达邻居（在距离阈值和时间阈值内）。接受参数index_center（整数: 给定核心点的索引）和df（dataframe: 单个车辆的轨迹点数据集）。函数返回给定核心点所有可达邻居点的索引集合。
fit(self, df)为实现ST-DBSCAN的聚类方法，用于接受某个车辆的轨迹点数据并完成STDBSCAN聚类。接受参数df（dataframe: 单个车辆的轨迹点数据集）。返回当前STDBSCAN类的实例本身。
class STDBSCAN(object): def __init__(self, spatial_threshold=500.0, temporal_threshold=30.0, min_neighbors=6): self.spatial_threshold = spatial_threshold self.temporal_threshold = temporal_threshold self.min_neighbors = min_neighbors self.labels_ = [] # 找到当前核心点的可达邻居 def retrieve_neighbors(self, index_center, df): neigborhood = [] # index_center为当前核心点索引，选取核心点对应的行数据 center_point = df.loc[index_center] # 根据时间阈值筛选可达点 min_time = center_point['collect_time'] - timedelta(minutes=self.temporal_threshold) max_time = center_point['collect_time'] + timedelta(minutes=self.temporal_threshold) df = df[(df['collect_time'] >= min_time) & (df['collect_time'] <= max_time)] # 根据距离阈值筛选可达点 for index, point in df.iterrows(): if index != index_center: distance = great_circle((center_point['lat'], center_point['lon']), (point['lat'], point['lon'])).m if distance <= self.spatial_threshold: neigborhood.append(index) # 返回所有可达点的索引 return neigborhood def fit(self, df): cluster_label = -1 noise = -1 unmarked = 777777 stack = [] # 丢弃原索引，保证新索引从0开始且连续 df = df.reset_index(drop=True) # 初始化聚类标签 self.labels_ = [unmarked] * df.shape[0] # 遍历所有轨迹点 for index, point in df.iterrows(): if self.labels_[index] == unmarked: neighborhood = self.retrieve_neighbors(index, df) if len(neighborhood) < self.min_neighbors: self.labels_[index] = noise else: # 发现核心点 # 为核心点分配新簇标签 cluster_label = cluster_label + 1 self.labels_[index] = cluster_label # 为该核心点的所有可达邻居点分配当前簇标签 for neighbor_index in neighborhood: self.labels_[neighbor_index] = cluster_label # 邻居点入栈 stack.append(neighbor_index) # 栈内点依次出栈 while len(stack) > 0: current_point_index = stack.pop() new_neighborhood = self.retrieve_neighbors(current_point_index, df) # 若current_point_index点也是核心点 if len(new_neighborhood) >= self.min_neighbors: for neighbor_index in new_neighborhood: neighbor_cluster = self.labels_[neighbor_index] if (neighbor_cluster != noise) & (neighbor_cluster == unmarked): self.labels_[neighbor_index] = cluster_label stack.append(neighbor_index) return self 3.进行聚类
在进行聚类之前，先使用shapely和geopy库实现了get_centermost_point函数。其输入数据cluster是列表类型，表示每一组聚类的点集。作用是在获得了每个聚类之后，计算出该聚类的中心点。
# 计算每个聚类的中心点 def get_centermost_point(cluster): # 计算整个点集合的质心点 centroid = (MultiPoint(cluster).centroid.x, MultiPoint(cluster).centroid.y) # 取点集合中离质心点最近的点为中心点 centermost_point = min(cluster, key=lambda point: great_circle(point, centroid).m) # 返回中心点 return tuple(centermost_point) 对单辆车的聚类函数cluster_traj，其输入数据data是dataframe类型，表示一辆车的轨迹数据。假设STDBSCAN类与cluster_traj函数放在同一文件中，对输入数据data进行一些预处理得到coords后，只需使用stdb = STDBSCAN(spatial_threshold=400.0, temporal_threshold=40.0, min_neighbors=4).fit(coords)就能完成聚类。怎么样，是不是很像scikit-learn中的DBSCAN用法？
# STDBSCAN聚类 def cluster_traj(data): # 创建副本 data = data.copy() # 提取dataframe中的经纬度列 coords = data[['collect_time', 'smoothed_lat', 'smoothed_lon']].copy() # 重命名列 coords.rename(columns={ 'smoothed_lat': 'lat', 'smoothed_lon': 'lon' }, inplace=True) # 将时间列转换为pandas的datetime格式用于后续计算 coords['collect_time'] = pd.to_datetime(coords['collect_time'].astype(str)) # 指定参数进行聚类，参数需根据数据集调整 stdb = STDBSCAN(spatial_threshold=400.0, temporal_threshold=40.0, min_neighbors=4).fit(coords) cluster_labels = stdb.labels_ # 将 cluster_labels 转换为 Pandas Series 并确保索引一致 cluster_labels = pd.Series(cluster_labels, index=coords.index) # 离群点的聚类标签为-1，其余数据聚成n类，标签为为0到n-1。num_clusters获得总共的聚类数n。 num_clusters = len(set(cluster_labels) - set([-1])) print('Clustered ' + str(len(data)) + ' points to ' + str(num_clusters) + ' clusters') # 将聚类标签添加为新列 c_label2 data['c_label2'] = cluster_labels # 聚类完成删除时间列并转换成数组，以便后续计算 coords = coords.drop(columns=['collect_time']).values # # 可选：输出聚类当前车的聚类情况 # clusters = pd.Series([coords[cluster_labels == n] for n in range(num_clusters)]) # print(clusters) # 输计算噪声点占总点数的比例 ratio = len(cluster_labels[cluster_labels[:] == -1]) / len(cluster_labels) print('噪声点占总点数的比例: ' + str(ratio)) # 只有聚类数量>1时才能计算指标 if num_clusters > 1: # 计算轮廓系数，作为聚类评价指标 sc_score = metrics.silhouette_score(coords, cluster_labels) print('轮廓系数: ' + str(sc_score)) # 计算DBI指标 dbi_score = metrics.davies_bouldin_score(coords, cluster_labels) print('戴维斯-布尔丁指数: ' + str(dbi_score)) print("\n") 4.单辆车的聚类评价
需要注意的是，每辆车聚类之后，还计算了噪声比、轮廓系数(SC)、戴维斯-布尔丁指数(DBI)来评价聚类效果。其中SC指标越接近1，聚类效果越好；DBI指标越小，聚类效果越好。
除此之外，还可以使用matplotlib库通过绘制散点图的方式，来肉眼观察这辆车的聚类效果，只需在cluster_traj函数中的return语句前插入下列代码（matplotlib绘制的散点图用于实验时判断聚类效果来调整参数，若要绘制更美观的图，可考虑使用folium库在地图上绘制轨迹点）：
# 获得每个聚类的中心点 centermost_points = clusters.map(get_centermost_point) # 将各个聚类的中心点存入rep_points lats, lons = zip(*centermost_points) rep_points = pd.DataFrame({'lon': lons, 'lat': lats}) # 绘制散点图 colors = list(mcolors.TABLEAU_COLORS.values()) # 使用Tableau颜色作为聚类颜色 noise_color = 'black' # 离群点颜色 fig, ax = plt.subplots(figsize=(12, 8)) for i, cluster in enumerate(clusters): if i == len(colors): # 如果聚类数超过颜色数，循环使用颜色 color = colors[i % len(colors)] else: color = colors[i] ax.scatter(cluster[:, 1], cluster[:, 0], s=30, c=color, marker='o', label='Cluster ' + str(i)) # 绘制离群点 noise_points = coords[cluster_labels == -1] ax.scatter(noise_points[:, 1], noise_points[:, 0], s=20, c=noise_color, marker='x', label='Noise points') ax.scatter(rep_points['lon'], rep_points['lat'], c='red', marker='*', s=100, label='Cluster Centers') ax.set_title('DBSCAN Clustering of Trajectory Data') ax.set_xlabel('Longitude') ax.set_ylabel('Latitude') ax.legend() plt.show() 二、整个数据集多辆车聚类
本人的数据集中包含多辆车的轨迹数据，这些数据统一存储一个CSV文件中，并且已经按照id和collect_time数据升序排序。
1.聚类
使用groupby的方式对车辆按id分组，每组分别调用cluster_traj即可。
#假设已经读入数据df clustered_data = pd.DataFrame() # 按车辆id分组，对每辆车的数据进行聚类 grouped = df.groupby('id') for name, group in grouped: print('车辆id：' + name + ' 轨迹点数：' + str(len(group))) clustered_group = cluster_traj(group) clustered_data = pd.concat([clustered_data, clustered_group], ignore_index=True) 2.整体评价
可以在函数外设置两个全局变量列表sc_scores和dbi_scores存储每辆车的评价指标。
# 全局变量用于存储指标 sc_scores = [] dbi_scores = [] 对cluser_traj函数中的计算轮廓系数部分添加sc_scores.append(sc_score)和dbi_scores.append(dbi_score)两行代码。即计算每辆车的评价指标的同时，将其加入外部的列表中。
if num_clusters > 1: # 计算轮廓系数，作为聚类评价指标 sc_score = metrics.silhouette_score(coords, cluster_labels) print('轮廓系数: ' + str(sc_score)) sc_scores.append(sc_score) # 计算DBI指标 dbi_score = metrics.davies_bouldin_score(coords, cluster_labels) print('戴维斯-布尔丁指数: ' + str(dbi_score)) dbi_scores.append(dbi_score) 通过sc_scores和dbi_scores两个列表的分析，例如求均值、中位数、画图查看分布等方式，可以评价整个数据聚类效果的好坏。
