基于DBACAN的道路轨迹点聚类

前言

很多针对道路轨迹的挖掘项目前期都需要对道路进行一段一段的分割成路段，然后对每一个路段来单独进行考察，如设定路段限速标识，超速概率等，如何对道路进行划分，其实是一个很有技巧性的活，最直白的有以下2种策略

• 道路栅格化

• 轨迹点聚类

下面分别对两种策略进行简单讲解。

道路栅格化

道路栅格化，简言之就是用一张纵横交错的网去尽可能覆盖道路所在的范围，这样，整个区域就被划分成一块一块的小矩形，形成栅格化，可以给每一个栅格编号，形成编号序列，而且可以判断出哪些栅格有轨迹点落入，哪些是没有轨迹点落入的，有轨迹点的栅格相对稀疏一些，此方法关键要考虑道路的经纬度最大范围和网眼大小，下面是道路栅格化处理主函数。

def roadRaster(road_data, unit_gap): #轨迹栅格化
    min_lng, max_lng = np.min(road_data['lng']), np.max(road_data['lng']) #经度范围
    min_lat, max_lat = np.min(road_data['lat']), np.max(road_data['lat']) #纬度范围
    lng_gap = max_lng - min_lng  
    lat_gap = max_lat - min_lat
    m = int(lng_gap/unit_gap)
    n = int(lat_gap/unit_gap)
    print(fleet_id,  min_lng, max_lng, min_lat, max_lat, m, n, (m-1)*(n-1))
    slice_lng = np.linspace(min_lng, max_lng, m)  #对经度等间距划分
    slice_lat = np.linspace(min_lat, max_lat, n) #对纬度等间距划分
    idx = 0
    for i in range(len(slice_lng)-1):
        for j in range(len(slice_lat)-1):
            raster_a_lng = slice_lng[i]
            raster_a_lat = slice_lat[j]
            raster_b_lng = slice_lng[i+1]
            raster_b_lat = slice_lat[j+1]
            idx +=1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

代码解读，首先，找出道路轨迹点经纬度最大最小值，然后对经纬度跨度进行等间距划分，然后对经纬度循环，不断生成栅格左下角点的经纬度对和右上角的经纬度对，由这样对顶角的点对就可以刻画出栅格，其中，unit_gap很关键，直接决定网眼大小，按下面经纬度小数点对应精度来粗略估计

轨迹聚类

轨迹聚类，就是根据历史行驶轨迹点的稠密程度来进行聚合成一簇一簇的轨迹点集合，其中同一簇的轨迹点尽可能靠在一起，不同一簇的轨迹点尽可能分散开来，然后把一簇的轨迹点范围提炼出来，如提取其四至，这样便把整个道路进行的切分。具体可以利用DBSCAN算法实现，DBACAN是一种基于密度的聚类算法，可以用于对道路轨迹点进行聚类。具体步骤如下：

• 初始化：将所有轨迹点标记为未访问状态，并设置一个固定的邻域半径r和最小聚类数量minPts。

• 随机选择一个未访问的点p，以p为中心，搜索其邻域内所有未访问点，并将这些点标记为已访问状态。

• 如果邻域内访问点的数量小于minPts，则将p标记为噪声点，否则创建一个新的聚类，并将p加入该聚类中。

• 遍历邻域内所有访问点的邻域，将其未访问的邻域点添加到聚类中，并将其标记为已访问状态。

• 重复2-4步，直到所有点都被访问过。

• 最后，将所有噪声点从聚类中去除。

需要注意的是，选择合适的邻域半径r和最小聚类数量minPts非常重要，这会影响到聚类结果的质量。可以通过试验不同的参数来获得最佳结果，下面是利用DBSCAN算法实现轨迹点聚类的主函数。

```python
def trajectoryCluster(trajectory,  k): # 使用KMeans， DBSCAN聚类算法进行路段划分
    ####先kmeans首次聚类，为每一个轨迹点分配一个类
    locations = np.array(trajectory[['lat','lng']]) #位置数据
    kmeans = KMeans(n_clusters = k)
    kmeans.fit(locations)
    labels = kmeans.labels_  #
    score = silhouette_score(locations, labels, metric='euclidean') #轮廓系数
    # print("一共聚了{}类, 轮廓系数为{}".format(labels.max() - labels.min(), score))
    cluster_label = pd.DataFrame({"cluster_label": labels})
    trajectory.reset_index(drop=True, inplace=True)
    cluster_label.reset_index(drop=True, inplace=True)
    cluster_data = pd.concat([trajectory, cluster_label], axis = 1, ignore_index=True) #带标签的行驶记录
    cluster_data.columns= ['lng', 'lat',  'cluster_label']
    cluster_data['cluster_label'] = [str(i) for i in cluster_data['cluster_label']]

    #####对每一个聚合出来的类，找出类代表，从而减少轨迹数量
    rep_trajectory = pd.pivot_table(cluster_data, index =['cluster_label'], values={'lng':'np.mean', 'lat':'np.mean'})
    rep_trajectory = rep_trajectory[['lng', 'lat']]
    # print("轨迹点代表\n", rep_trajectory)
    return rep_trajectory

def calDistance(point1, point2): #计算两点之间的曼哈顿距离
    manhattan_distance = np.abs(point1[0] - point2[0]) + np.abs(point1[1]-point2[1])
    return manhattan_distance

def trajectorySort(trajectory, init_point): #对无序的轨迹点进行排序得出路径
    # print("初始位置", trajectory[0])
    num_points = len(trajectory) #总轨迹点数
    visited = [False]*num_points
    path = [] #用来存放整理后的轨迹点
    
    current_point = init_point #trajectory[0]  #当前轨迹点
    path.append(current_point) #把当前轨迹点追加近路径
    visited[0] = True

    while len(path)< num_points:
        min_distance = float('inf')
        nearest_point = None

        for i, point in enumerate(trajectory):
            if  not visited[i]:
                distance = calDistance(current_point, point)
                if distance < min_distance:
                    min_distance = distance
                    nearest_point = point

        if nearest_point is not None:
            path.append(nearest_point)
            visited[trajectory.index(nearest_point)] = True
            current_point = nearest_point
    # path.append(path[0]) #形成回路
    return path

################重要分割线#######################

def myScore(estimator, X): #轮廓系数
    labels = estimator.fit_predict(X)
    score = silhouette_score(X, labels, metric='euclidean')
    return score

def roadCluster(trajectory): # 使用DBSCAN聚类算法进行路段划分
    sample_num = int(0.6*len(trajectory))
    print(sample_num)
    trajectory_sample =  trajectory.sample(sample_num) #随机抽样60%样本点 
    locations = np.array(trajectory_sample[['lat','lng']]) #位置数据
    param_grid = {"eps":[0.0005, 0.001, 0.003,  0.005, 0.006, 0.01],
                  "min_samples":[6,  9, 15, 20, 30, 50, 70]
                  } # epsilon控制聚类的距离阈值，min_samples控制形成簇的最小样本数
    dbscan = DBSCAN()
    grid_search = GridSearchCV(estimator= dbscan, param_grid=param_grid, scoring=myScore)
    grid_search.fit(locations)
    print("best parameters:{}".format(grid_search.best_params_))
    print("label:{}".format(grid_search.best_estimator_.labels_))
    labels = grid_search.best_estimator_.labels_  #-1表示离群点
    score = silhouette_score(locations, labels, metric='euclidean') #轮廓系数
    total_cluster = labels.max() - labels.min()
    print("一共聚了{}类, 轮廓系数为{}".format(total_cluster, score))
    road_label = pd.DataFrame({"road_label": labels})
    trajectory_sample.reset_index(drop=True, inplace=True)
    road_label.reset_index(drop=True, inplace=True)
    cluster_data = pd.concat([trajectory_sample, road_label], axis = 1, ignore_index=True) #带标签的行驶记录
    cluster_data.columns= ['lng', 'lat', 'speed', 'road_label']
    cluster_data['road_label'] = [str(i) for i in cluster_data['road_label']]
    print(cluster_data)
    return cluster_data
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

代码解读，聚类的对象locations由经纬度组成的2维数组，通过grid_search 来寻找最佳的超参数epsilon和min_samples,最后把标签类和原先的轨迹拼接起来，相当于给原先的每一个轨迹点打一个类别标签，聚类后，可以只提炼出每一类的经纬度中位数进行可视化，即用一个点来代表这一簇，效果会显得更加稀疏明显

参考资料

1，经纬度坐标小数位与精度的对应关系
https://blog.csdn.net/lang_niu/article/details/123550453

2，基于DBSCAN算法的营运车辆超速点聚类分析
https://max.book118.com/html/2018/0407/160435287.shtm