2024年华中杯数学建模B题使用行车轨迹估计交通信号灯周期问题解题全过程文档及程序

2024年华中杯数学建模

B题使用行车轨迹估计交通信号灯周期问题

原题再现

某电子地图服务商希望获取城市路网中所有交通信号灯的红绿周期，以便为司机提供更好的导航服务。由于许多信号灯未接入网络，无法直接从交通管理部门获取所有信号灯的数据，也不可能在所有路口安排人工读取信号灯周期信息。所以，该公司计划使用大量客户的行车轨迹数据估计交通信号灯的周期。请帮助该公司解决这一问题，完成以下任务。已知所有信号灯只有红、绿两种状态。
1. 若信号灯周期固定不变，且已知所有车辆的行车轨迹，建立模型，利用车辆行车轨迹数据估计信号灯的红绿周期。附件1中是5个不相关路口各自一个方向连续1小时内车辆的轨迹数据，尝试求出这些路口相应方向的信号灯周期，并按格式要求填入表1。
2. 实际上，只有部分用户使用该公司的产品，即只能获取部分样本车辆的行车轨迹。同时，受各种因素的影响，轨迹数据存在定位误差，误差大小未知。讨论样本车辆比例、车流量、定位误差等因素对上述模型估计精度的影响。附件2中是另外5个不相关路口各自一个方向连续1小时内样本车辆的轨迹数据，尝试求出这些路口相应方向的信号灯周期，按同样的格式要求填入表2。
3. 如果信号灯周期有可能发生变化，能否尽快检测出这种变化，以及变化后的新周期？附件3中是另外6个不相关路口各自一个方向连续2小时内样本车辆的轨迹数据，判断这些路口相应方向的信号灯周期在这段时间内是否有变化，尝试求出周期切换的时刻，以及新旧周期参数，按格式要求填入表3，并指明识别出周期变化所需的时间和条件。
4. 附件4是某路口连续2小时内所有方向样本车辆的轨迹数据，请尝试
识别出该路口信号灯的周期。
 附件1：路口A1、A2、A3、A4、A5各自一个方向连续1小时内车辆轨迹数据
 附件2：路口B1、B2、B3、B4、B5各自一个方向连续1小时内样本车辆轨迹数据
 附件3：路口C1、C2、C3、C4、C5、C6各自一个方向连续2小时内样本车辆轨迹数据
 附件4：路口D所有方向连续2小时内样本车辆轨迹数据
 附件5：数据文件说明及结果表格
1、轨迹数据文件格式。适用于附件1-附件4所有轨迹数据文件。纯文本文件，第一行为标题行，各列以英文逗号分隔，共5列，分别为时间点、车辆ID、当前位置X坐标、当前位置Y坐标。时间点单位为秒，第0秒开始，每1秒采样一次。坐标单位为米。车辆ID仅用于区分同一个文件中的不同车辆。车辆ID不一定是连续编号。不同文件中，相同ID的车辆没有任何联系。同一车道可能只允许一个方向前进，也可能允许两个方向前进，如直行或左转、直行或右转等。
2、表1：路口A1-A5各自一个方向信号灯周期识别结果
在这里插入图片描述
3、表2：路口B1-B5各自一个方向信号灯周期识别结果

4、表3：路口C1-C6各自一个方向信号灯周期识别结果

说明：“周期切换时刻”是指信号灯周期发生变化的具体时间点，以第一个变化后的时长区间的起点计。如果信号灯周期没有变化，则“周期切换时刻”填写“无”。如果信号灯周期多次切换，按照上述格式，自行延长表格依次填写。

整体求解过程概述(摘要)

本文根据某电子地图服务商提供的行车轨迹数据，建立了基于两阶段STFT的时频分析模型和基于DBSCAN聚类分析改进的周期检测模型对不同路口的周期进行估计；然后利用贝叶斯变点检测求出周期变化时各阶段的周期参数。
问题一中，本文首先对 A1-A5 五个互不相关路口的行车轨迹数据进行离群点检测、轨迹连续性检验，然后进行异常值处理。然后，通过指标构建得到车辆运行加速度状态的时序数据，本文建立了基于两阶段STFT的时频分析模型，对时序数据进行周期特征挖掘，得出信号灯的周期性规律。然后通过绘制功率谱密度图得到信号灯完整周期长度，再将车辆处于停车状态最长的时间作为红灯周期，进而得到绿灯周期。最终得到 A1-A5 路口的红灯周期分别为[73,59,81,69,63];绿灯周期分别为[33,28,26,19,26]。
问题二中，由于样本车辆减少，存在定位偏差等情况，因此，本文在第一问所建立模型的基础上融入DBSCAN聚类分析，得到红绿灯完整周期时长或者完整周期时长的倍数，然后利用基于两阶段STFT的时频分析模型对红绿灯周期进行估计，得到B1-B5 路口的红灯周期分别为[78,83,54,78,95]，绿灯周期分别为[27,33,34,27,21]。最后对模型设计对比实验分析探究影响因素，发现车流量减小对模型的影响较大，定位误差因素对模型基本无影响，发现模型有较高的鲁棒性。
问题三中，每个路口信号灯周期均可能发生变化，故本文使用贝叶斯变点检测对周期的变化点进行检验；然后将时序数据沿变化点进行分段，对分段后的每段时序数据使用加Manning窗的FFT算法求解周期，得到C1-C6路口的主要完整周期长度分别为[88,88,105,105,88,105]，每个路口周期变化次数分别为[6,5,3,0,4,0]。
问题四中，附件4给出的是某路口D所有方向的车辆轨迹数据，本文首先建立多方向汽车轨迹分解模型对该路口的车辆行驶方向进行分解，得到沿一个方向的所有车辆轨迹；再使用前文模型求解各个方向的周期，最后使用滑动窗口稳态检测对周期进行检验，得到路口所以方向完整红绿灯周期时长均为142s，且对向车道信号灯周期时间相同，符合现实红绿灯规律。

模型假设：

1、假设每个路口在正常情况下,不会发生交通事故、堵车等意外情况,道路通行仅受红绿灯信号的影响；
2、针对附件提供的数据,假设这些数据是真实可靠的,并且车辆在行驶过程中均遵守交通规则,因此可以直接基于这些数据进行计算分析；
3、假设题目中给出的若干路口的红绿灯信号情况是相互独立的,即上一个路口的红绿灯状态不会对下一个路口的红绿灯信号产生任何影响。

问题分析：

针对问题一，首先需要对车辆行车轨迹数据进行清洗，剔除异常值。这可以通过确实轨迹正方向、离群点检验、轨迹连续性检验等方式来识别异常点，并用通过差值进行替换，确保数据的连贯性。接下来,构建一系列反映车辆运行状态的指标,包括车辆运行状态、车辆速度、车辆加速度等。这些指标能够共同形成车辆状态的时序数据。最后,利用基于两阶段STFT的时频分析模型，可以对这些时序数据进行周期特征挖掘,得出车辆行驶的周期性规律。

针对问题二,只有部分用户在使用该公司的产品，因此只能获取部分样本车辆的行车轨迹数据。而且,样本车辆比例的减少、车流量的减少以及存在定位误差等因素，都会对上述STFT预测模型的预测精度产生较大的影响。因此，本文在第一问的基础上结合DBSCAN聚类分析，通过聚类可以得到这一次红灯结束的时间到下一次红灯结束的时间，即为红绿灯完整周期时长或者周期时长的倍数，进而得到红绿灯完整周期得大致范围，从而对第一问中的两阶段STFT的时频分析模型进行修正。

针对问题三，贝叶斯变点检测是一种用于发现时序数据中周期性变化点的有效方法。它建立在概率模型的基础之上,利用贝叶斯推断来识别数据序列中的变化点。具体而言,贝叶斯变点检测假设数据序列可以由不同的概率模型描述,每个模型对应于序列的一个状态。通过计算不同模型的后验概率,可以确定序列中最可能出现变化点的位置。与其他变点检测方法相比,贝叶斯变点检测不仅能够检测出变化的时间点,还可以量化变化的概率,从而更好地反映变化的程度和确定性。因此本文使用贝叶斯变点检测检测出周期变化点，然后将时序数据沿变化点进行分段，对每段时序数据使用加manning窗的FFT算法求解周期，通过上述方式得到周期的变化时间点，然后对时间点分割，对分割后的数据采用前文模型，得到具体周期时间。

针对问题四，问题四涉及某路口的全方位数据，我们可以将其拆分为前述问题的子问题。具体来说，先对该路口的车辆运行模式进行分解分析，然后应用前文提及的模型对各个方向进行周期性特征提取。接下来，可以利用改进的线性互相关算法对周期预估结果进行验证。通过这种方式，们可以将问题四转化为先前提出的问题类型，并继续采用前文介绍的分析方法进行求解。这样不仅避免重复阐述，也能充分利用前期的建模和分析思路。

模型的建立与求解整体论文缩略图

在这里插入图片描述

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

部分程序代码：

import pandas as pd 
import matplotlib.pyplot as plt 
from scipy import stats 
import numpy as np 
filename = 'data/2/B5.csv' 
df = pd.read_csv(filename) 
#对df.x,df.y做散点图，并做离群点检测 
plt.figure(figsize=(8, 6)) 
plt.rcParams['font.sans-serif'] = ['SimHei'] 
plt.rcParams['axes.unicode_minus'] = False 
#设置字体大小 
plt.rcParams['font.size'] = 14 
plt.scatter(df['x'], df['y'], label='车辆轨迹数据') 
plt.xlabel('x',fontsize = 16) 
plt.ylabel('y',fontsize = 16) 
plt.title('车辆轨迹散点图') 
###只放一部分
    return outlier_indices 
 
outlier_indices = detect_outliers(df['x'], df['y']) 
outliers = df.iloc[outlier_indices]
print(f"Outlier indices: {outlier_indices}") 
print(f"Number of outliers: {len(outlier_indices)}") 
print(outliers) 
#在图中画出异常点 
plt.figure(figsize=(8, 6)) 
plt.scatter(df['x'], df['y'], label='车辆轨迹数据') 
plt.scatter(outliers['x'], outliers['y'], color='red', label='异常点') 
plt.xlabel('x') 
plt.ylabel('y') 
plt.title('车辆轨迹散点图') 
plt.legend() 
plt.savefig('data2/img/B5_scatterDetect.png') 
plt.show() 
 
 
from statsmodels.tsa.seasonal import seasonal_decompose 
# 
#读取data中1中的A1.csv文件 
import pandas as pd 
import math 
import os 
import matplotlib.pyplot as plt 
import seaborn as sns 
import numpy as np 
import pandas as pd 
 2 
    df.loc[i[1].index,'x_y'] = df.loc[i[1].index,'x_y'].fillna(2) 
#对df中x_y列中的值大于0.5的行，添加一个标签，值为1，否则为0，列名为label 
) 
ts = data['count'] 
#探究ts的周期性 
decomposition = seasonal_decompose(ts,period=72, 
                                model="add") 
trend = decomposition.trend 
seasonal = decomposition.seasonal 
residual = decomposition.resid 
plt.plot(ts, label='Original') 
plt.legend(loc='best') 
plt.plot(trend, label='Trend') 
plt.legend(loc='best') 
plt.savefig('data2/img/D_trend.png',dpi=200) 
plt.show() 
 
 
import numpy as np 
import pandas as pd 
import matplotlib.pyplot as plt 
from statsmodels.graphics.tsaplots import plot_acf 
 
plt.legend() 
plt.show() 
 
# 自相关函数（ACF）分析 
plt.figure(figsize=(12, 4)) 
plot_acf(diff_data, lags=1, title="差分序列自相关函数（ACF）") 
plt.show() 
 
 
array = [ 134,  246,  398,  485,  685,  948, 1028, 1213, 1653, 1828, 2004, 
        2092, 2245, 2621, 2796, 2886, 3148, 3325, 3390, 3500, 3589, 3677, 
        3764, 3940, 4028, 4292, 4380, 4469, 4646, 4820, 4997, 5172, 5436, 
        5612, 5704, 5964, 6140, 6228, 6381, 6581, 6828, 7108, 7196] 
 
percentage_threshold = 8  # 百分比阈值 
pe = 88 
differences = [] 
previous_difference = array[1] - array[0] 
change_points = [] 
 
for i in range(2, len(array)): 
    current_difference = array[i] - array[i - 1] 
    r
    sample_freq = fftfreq(fft_series.size) 
    pos_mask = np.where(sample_freq > 0) 
    freqs = sample_freq[pos_mask] 
    powers = power[pos_mask] 
    top_k_seasons = 3 
    top_k_idxs = np.argpartition(powers, -top_k_seasons)[-top_k_seasons:] 
    top_k_power = powers[top_k_idxs] 
    fft_periods = (1 / freqs[top_k_idxs]).astype(int) 
 
    print(f"top_k_power: {top_k_power}") 
    print(f"fft_periods: {fft_periods}") 
 
    NFFT = 256 
    TEMP = np.abs(np.fft.fft(ts.values, n=NFFT)) 
    PSD = (TEMP[0:NFFT//2+1])**2 / NFFT 
    frequencies = np.linspace(0, 0.78/2, NFFT//2+1)  # 单边频率轴，从0到采样率一半 
    periods = 1./frequencies  # 将频率转换为周期 
 fig, ax = plt.subplots() 
    ax.semilogy(periods, PSD, linewidth=1.5)  # 使用对数尺度展示功率谱密度 
    #标出periods范围为50~100间的对应PSD的最大值点 
    max_index = np.argmax((periods >= 50) & (periods <= 100) & (PSD > 0))  # 需确保 PSD 大于 
0，避免选取到噪声 
    ax.scatter(periods[max_index], PSD[max_index], color='red', marker='o', label=f"最大PSD点对
应周期为: {periods[max_index]:.2f}") 
    # print(periods[np.argmax(PSD)], PSD[np.argmax(PSD)]) 
    # ax.scatter(periods[np.argmax(PSD)], PSD[np.argmax(PSD)], color='red', s=100) 
    ax.set_xlabel('周期（秒）') 
    ax.set_ylabel('PSD值') 
    ax.grid(True) 
    ax.legend() 
    plt.rcParams['font.sans-serif'] = ['SimHei'] 
    plt.rcParams['axes.unicode_minus'] = False 
    plt.savefig('data2/img/STFT.png',dpi=200) 
    plt.show() 
 
    from statsmodels.tsa.stattools import acf 
 
    for lag in fft_periods: 
        acf_score = acf(ts.values, nlags=lag)[-1] 
        print(f"lag: {lag} fft acf: {acf_score}") 
 
    expected_lags = np.array([30, 60, 90]).astype(int) 
    for lag in expected_lags: 
        acf_score = acf(ts.values, nlags=lag, fft=False)[-1] 
        print(f"lag: {lag} expected acf: {acf_score}") 
    pl = df 
    print(pl.groupby('vehicle_id')['state'].sum().max(),pl.groupby('vehicle_id')['state'].sum().idxmax()) 
    print(pl.groupby('vehicle_id')['state'].sum().nlargest(2).iloc[1],pl.groupby('vehicle_id')['state'].sum
 ().nlargest(2).index[1]) 
    print(pl.groupby('vehicle_id')['state'].sum().nlargest(3).mean()) 
    sns.distplot(pl.groupby('vehicle_id')['state'].sum(),bins=20) 
STFT('data/4/D.csv')