Python 3.12 + Pandas 复现 MathorCup B题：共享单车时空分布与OD矩阵计算

2026/7/8 23:10:03

Python 3.12 + Pandas 复现 MathorCup B题：共享单车时空分布与OD矩阵计算实战指南

1. 项目背景与数据准备

共享单车作为城市短途出行的重要工具，其运营效率直接影响用户体验和企业效益。2017年MathorCup数学建模B题以共享单车调度优化为背景，要求参赛者分析时空分布特征并构建OD矩阵。本文将使用Python 3.12和Pandas生态系统完整复现该案例。

核心数据集应包含以下字段：

start_time: 骑行开始时间（datetime格式）
end_time: 骑行结束时间（datetime格式）
start_region: 出发区域编码（整数）
end_region: 到达区域编码（整数）
user_id: 用户唯一标识符

典型数据加载与预处理代码：

import pandas as pd import numpy as np # 加载原始数据 raw_data = pd.read_csv('bike_trips.csv', parse_dates=['start_time', 'end_time'], dtype={'start_region': 'int8', 'end_region': 'int8'}) # 计算骑行时长（分钟） raw_data['duration'] = (raw_data['end_time'] - raw_data['start_time']).dt.total_seconds() / 60 # 过滤异常数据 clean_data = raw_data[ (raw_data['duration'] > 1) & (raw_data['duration'] < 180) & (raw_data['start_region'].between(1, 10)) & (raw_data['end_region'].between(1, 10)) ].copy()

提示：实际项目中建议使用pandas-profiling或ydata-profiling生成数据质量报告，快速识别异常值和缺失值。

2. 时空分布特征分析

2.1 时间维度分析

小时级骑行量统计：

# 提取小时特征 clean_data['hour'] = clean_data['start_time'].dt.hour # 分时段统计 hourly_stats = clean_data.groupby('hour').agg( trip_count=('user_id', 'count'), avg_duration=('duration', 'mean') ).reset_index()

可视化时段分布（需配合Matplotlib/Seaborn）：

import matplotlib.pyplot as plt import seaborn as sns plt.figure(figsize=(12, 6)) sns.lineplot(data=hourly_stats, x='hour', y='trip_count', marker='o') plt.title('Hourly Bike Trip Distribution') plt.xlabel('Hour of Day') plt.ylabel('Trip Count') plt.grid(True) plt.show()

典型时间分布特征：

早高峰：7:00-9:00（通勤需求）
午间小高峰：11:00-13:00（午餐出行）
晚高峰：17:00-19:00（下班通勤）
夜间低谷：22:00-次日5:00

2.2 空间维度分析

区域热度排名：

departure_rank = clean_data['start_region'].value_counts().sort_index() arrival_rank = clean_data['end_region'].value_counts().sort_index() region_stats = pd.DataFrame({ 'departures': departure_rank, 'arrivals': arrival_rank, 'net_flow': arrival_rank - departure_rank })

区域类型	特征描述
净流入区	商业中心、交通枢纽
净流出区	住宅区、学校
平衡区	混合功能区

3. OD矩阵构建与优化

3.1 基础OD矩阵计算

# 构建原始OD矩阵 od_matrix = clean_data.groupby(['start_region', 'end_region']).size().unstack(fill_value=0) # 标准化处理（出行比例） od_matrix_normalized = od_matrix.div(od_matrix.sum(axis=1), axis=0)

3.2 时空增强OD矩阵

分时段OD矩阵示例：

def build_time_slice_od(df, time_ranges): results = [] for name, (start_h, end_h) in time_ranges.items(): slice_data = df[df['hour'].between(start_h, end_h)] od = slice_data.groupby(['start_region', 'end_region']).size().unstack(fill_value=0) od['time_slice'] = name results.append(od) return pd.concat(results) time_ranges = { 'morning': (7, 9), 'noon': (11, 13), 'evening': (17, 19) } enhanced_od = build_time_slice_od(clean_data, time_ranges)

3.3 动态OD矩阵优化

基于滑动窗口的实时更新：

def dynamic_od_update(window_minutes=30): time_range = pd.date_range( start=clean_data['start_time'].min().floor('H'), end=clean_data['start_time'].max().ceil('H'), freq=f'{window_minutes}T' ) dynamic_ods = [] for i in range(len(time_range)-1): window_data = clean_data[ (clean_data['start_time'] >= time_range[i]) & (clean_data['start_time'] < time_range[i+1]) ] od = window_data.groupby(['start_region', 'end_region']).size().unstack(fill_value=0) od['window_start'] = time_range[i] dynamic_ods.append(od) return pd.concat(dynamic_ods)

4. 调度优化模型实现

4.1 需求预测模型

from sklearn.ensemble import RandomForestRegressor from sklearn.model_selection import train_test_split # 特征工程 features = pd.get_dummies(clean_data[['hour', 'start_region', 'end_region']]) target = clean_data['duration'] # 训练测试分割 X_train, X_test, y_train, y_test = train_test_split( features, target, test_size=0.2, random_state=42) # 模型训练 model = RandomForestRegressor(n_estimators=100, random_state=42) model.fit(X_train, y_train)

4.2 调度路径优化

使用OR-Tools实现VRP：

from ortools.constraint_solver import routing_enums_pb2 from ortools.constraint_solver import pywrapcp def optimize_route(demand_matrix, distance_matrix): # 创建距离回调函数 def distance_callback(from_index, to_index): from_node = manager.IndexToNode(from_index) to_node = manager.IndexToNode(to_index) return distance_matrix[from_node][to_node] # 初始化路由模型 manager = pywrapcp.RoutingIndexManager( len(distance_matrix), 1, 0) # 1辆车，0为仓库 routing = pywrapcp.RoutingModel(manager) transit_callback_index = routing.RegisterTransitCallback(distance_callback) routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index) # 设置搜索参数 search_parameters = pywrapcp.DefaultRoutingSearchParameters() search_parameters.first_solution_strategy = ( routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC) # 求解 solution = routing.SolveWithParameters(search_parameters) # 提取最优路径 index = routing.Start(0) route = [] while not routing.IsEnd(index): route.append(manager.IndexToNode(index)) index = solution.Value(routing.NextVar(index)) route.append(manager.IndexToNode(index)) return route

5. 可视化与结果分析

5.1 热力图可视化

import plotly.express as px fig = px.imshow(od_matrix_normalized, labels=dict(x="Destination", y="Origin", color="Trip Ratio"), x=[f'Zone {i}' for i in range(1, 11)], y=[f'Zone {i}' for i in range(1, 11)], color_continuous_scale='Viridis') fig.update_layout(title='Normalized OD Matrix') fig.show()

5.2 动态调度模拟

def simulate_dynamic_scheduling(od_series, initial_inventory): inventory = initial_inventory.copy() records = [] for ts, od in od_series.items(): # 记录当前库存 record = inventory.copy() record['timestamp'] = ts records.append(record) # 模拟骑行消耗 for origin in od.index: for dest in od.columns: trips = od.loc[origin, dest] inventory[origin] -= trips inventory[dest] += trips # 库存修正（防止负值） inventory = inventory.clip(lower=0) return pd.DataFrame(records).set_index('timestamp')

6. 工程化实践建议

性能优化技巧：
- 对大规模数据使用pandas.DataFrame.eval()进行链式运算
- 使用numba加速数值计算密集型函数
- 对分类数据使用category数据类型减少内存占用

部署架构：

graph TD A[实时数据流] --> B(Kafka消息队列) B --> C{Spark流处理} C --> D[OD矩阵计算] C --> E[需求预测] D --> F[调度优化引擎] E --> F F --> G[调度指令下发]

异常处理模式：

class BikeDataException(Exception): pass def validate_od_matrix(matrix): if matrix.isnull().any().any(): raise BikeDataException("OD矩阵包含空值") if (matrix < 0).any().any(): raise BikeDataException("OD矩阵包含负值") return True

在实际项目中，我们还需要考虑：