大数据正在渗透各行各业，甚至能跟你考试能力测试、患上某种疾病的机率等非常生活化的场景应用都发生紧密的联系。今后大数据在我们的生活中就像是水和电一样，让社会整个信息质量更好、让信息利用效率更高效。本文揭示了新加坡政府是如何利用大数据技术来捕获引发地铁被中断的反常列车，我们得以再一次见识大数据技术的神奇力量。
世界著名未来学家托夫勒曾说改变这个世界的力量有三种暴力、知识、金钱，而如今我们的世界正在被第四种力量改变，那就是大数据!大数据不管应用在哪个行业它的核心都是通过技术来获知事情发展的真相，最终利用这个“真相”来更加合理的配置资源。具体来说，要实现大数据的核心价值，还需要前两个重要的步骤，第一步是通过“众包”的形式收集海量数据，第二步是通过大数据的技术途径进行“全量数据挖掘”，最后利用分析结果进行“资源优化配置”。说白了，大数据最终的落地就是资源优化配置。
本文揭示了新加坡政府是如何利用大数据技术来捕获引发地铁被中断的反常列车，我们得以再一次见识大数据技术的神奇力量。
最近几个月，新加坡的地铁环线(MRT Circle Line)遭到了一连串的神秘中断，对数以千计的乘客造成了很大的混乱和痛苦。
同大多数同事一样，我每天早晨搭乘环线地铁到办公室。因此，在11月5日，当我所在的团队有调查原因的机会时，我就毫不犹豫地自告奋勇参加了。
根据新加坡地铁公司(SMRT)和新加坡陆路交通管理局(Land Transport Authority，LTA)的先前调查，我们知道这些事件是由于某种形式的信号干扰造成的，导致了一些列车的信号丢失。信号丢失会触发那些列车中的制动安全功能，并使它们沿着轨道随机停止。
但是这起第一次发生在八月份的事件——似乎是随机发生的，使调查小组很难找到确切的原因。
我们获得了由SMRT编译的数据集，其中包含以下信息：
- 每个事件的日期和时间
- 事件的位置
- 涉及的列车编号
- 列车的方向
我们开始清理数据，在Jupyter Notebook中进行工作，这是一个流行的编写和记录Python代码的工具。
像往常一样，第一步是导入一些有用的Python库。
[mw_shl_code=python,true]import math
import xlrd
import itertools as it
import numpy as np
import pandas as pd
from datetime import datetime [/mw_shl_code]
片段1
然后我们从原始数据中提取有用的部分。
[mw_shl_code=python,true]dfincidents_0 = pd.read_excel('CCL EVAC E-brake occurrences hourly update_mod.xlsx', sheetname='Aug Sep')
dfincidents_1 = pd.read_excel('CCL EVAC E-brake occurrences hourly update_mod.xlsx', sheetname='Nov')
# Incident data for Nov had different format
dfincidents_1['Time'] = dfincidents_1['Time'].str.strip('hrs').str.strip(' ')
dfincidents_1['Time'] = pd.to_datetime(dfincidents_1['Time'], format='%H%M').dt.time
dfincidents = pd.concat([dfincidents_0, dfincidents_1])
# Reset the index because they were concatenated from two data sources
dfincidents.reset_index(inplace=True, drop=True) [/mw_shl_code]
片段2
我们将日期和时间列合并为一个标准列，以便更容易地将数据可视化：
[mw_shl_code=python,true]def datetime_from_date_and_time(row):
"""
Combines the date column and time column into a single column
"""
d = row['Date']
t = row['Time']
return datetime(
d.year, d.month, d.day,
t.hour, t.minute, t.second
)
# Add DateTime to the data for easier visualization
dfincidents['DateTime'] = dfincidents.apply(datetime_from_date_and_time, axis=1) [/mw_shl_code]
片段3
这就产生了如下表格：
截图1：初始处理的输出
最初的可视化没有明确的答案
我们在初步探索性的分析中找不到任何明显的答案，如下图所示：
1.事件发生在一天之中，整天的事件数量反映了高峰和非高峰旅行时间。
图1 出现反映高峰和非高峰旅行时间的次数
2.事故发生在环线上的各个地点，西侧发生的事件略多一些。
图2 干扰的原因似乎与位置无关
3.只有一辆或两辆列车的话，信号干扰并不会产生影响，但在这条环线上有许多列车。“PV”是“客运车辆”(Passenger Vehicle)的缩写。
图3 60辆不同列车受到信号干扰
Marey图表：显示时间、位置和方向
我们的下一步是将多维度纳入探索性分析。
我们的灵感来自Marey图表，这是1983年Edward Tufte出版的经典著作《定量信息的视觉显示》(《The Visual Display of Quantitative Information》)最近，它被Mike Barry和Brian Card应用在波士顿地铁系统的大规模可视化项目：
截图2 摘自http://mbtaviz.github.io/
在该图表中，垂直轴表示时间——按时间顺序从上到下;而水平轴表示沿着列车线路的车站。对角线则表示列车运行。
在我们的Marey图表中开始绘制轴：
图4 环线版本的一个空白Marey图表
在正常情况下，在HarbourFront和Dhoby Ghaut之间运行的列车将在类似与此的线路上移动，每次单程行程只需要一个小时以上：
图5 环线上列车运行的程式化表示
我们的目的是在这个图表上绘制事件——是点而不是线。
准备可视化数据
首先，我们将站名从三字母代码转换为数字：
- Marina Bay到Promenade之前：0到1.5;
- Dhoby Ghaut到HarbourFront：2到29。
如果事故发生在两个站之间，则将其表示为0.5加上两个站号中的较小数。例如，如果事件发生在HarbourFront(29)和Telok Blangah(28)之间，则位置将是“28.5”。这样我们就很容易绘制沿水平轴的点。
[mw_shl_code=python,true]stations=("MRB,BFT,DBG,BBS,EPN,PMN,NCH,SDM,MBT,DKT,PYL,MPS,TSG,BLY,SER,"
"LRC,BSH,MRM,CDT,BTN,FRR,HLV,BNV,ONH,KRG,HPV,PPJ,LBD,TLB,HBF").split(',')
def loc_id(station1, station2 = None):
"""
Translates a 3-letter station code to a number,
or a pair of 3-letter station codes to a number.
Single stations are represented as whole numbers.
Locations between stations are represented with a .5.
Example:
loc_id('MRB') # 0 (Marina Bay)
loc_id('MRB', 'BFT') # 0.5 (Between Marina Bay and Bayfront)
loc_id('DBG') # 2 (Dhoby Ghaut)
loc_id('HBF') # 29
loc_id('HBF', 'TLB') # 28.5 (Between Harbourfront and Telok Blangah)
loc_id('HBF', 'DBG') # throws and error, because these stations are not adjacent
"""
if station2 == None or station2 == 'nan' or (type(station2) is float and math.isnan(station2)):
# Single stations
return stations.index(station1)
else: # Pairs of stations -- take the average to get the 0.5
stn1_index = stations.index(station1)
stn2_index = stations.index(station2)
# Handle the branch at Promenade
if (set(['PMN', 'EPN']) == set([station1, station2])):
return float(stations.index('EPN')) + 0.5
elif set(['PMN', 'BFT']) == set([station1, station2]):
return float(stations.index('BFT')) + 0.5
else:
# Require station pairs to be adjacent stations
assert(math.fabs(stn1_index - stn2_index) == 1)
return float(stn1_index + stn2_index) / 2 [/mw_shl_code]
片段 4
然后我们计算位置ID的数字……
[mw_shl_code=python,true]def loc_id_from_stations(row):
try:
# This handles entries with both "Station from" and "Station to"
# and entries with only "Station from"
return loc_id(row['Station from'], str(row['Station to']))
except ValueError:
# Some entries only have "Station to" but no
# "Station from"
return loc_id(row['Station to']) [/mw_shl_code]
片段 5
并添加到数据集：
[mw_shl_code=python,true]# Select only some columns that we are interested in
sel_dfincidents = dfincidents[['DateTime', 'PV', 'Bound', 'Station from', 'Station to', 'Event', 'Remarks']]
# Add the location ID into the dataset
sel_dfincidents['LocID'] = sel_dfincidents.apply(loc_id_from_stations, axis=1) [/mw_shl_code]
片段6
然后我们得到如下表格：
截图3 添加位置ID后的输出表
通过数据处理，我们能够创建所有紧急制动事件的散点图。这里的每个点代表一个事件。我们还是无法发现任何明显的事故模式。
图6 信号干扰事件表示为散点图
接下来，我们通过将每个事件表示为指向左侧或右侧的三角形，而不是点，将列车方向添加到图表中：
图7 方向由箭头和颜色表示。
它看起来相当随机，但当我们放大到如下的图表，一个模式似乎似乎浮出了水面：
图8 上午6点到10点之间的事件
如果你仔细阅读图表，你会注意到，故障似乎按顺序发生。当一趟列车受到干扰时，另一趟在同一方向行驶的列车很快就会受到波及。
信号如何互相干扰?
在这一点上依然不明，一趟列车是罪魁祸首。
我们已经证实的是，似乎有一个随时间和位置相关的模式：事件一个接一个地发生，与上一个事件的方向相反。似乎有一条“破坏的痕迹”，它会不会是导致数据集中那些事件的诱因?
事实上，连接事件的假想线看上去与截图2的Marey图表可疑地类似。干扰的原因会不会是对面轨道的列车?
图9 它会是一趟相反方向行驶的列车吗?
我们决定检验这个“反常列车”的假说。
我们已经知道，沿着环线的车站之间的行驶时间在两到四分钟之间。这意味着如果发生4分钟的间隔，我们可以将所有紧急制动事件分组在一起。
[mw_shl_code=python,true]def same_cascade(i, j):
"""
Given a pair of incidents (i,j), returns true if:
t <= d * 4 mins
where t is the time difference between occurrences
and d is the distance (measured by difference in location ID).
Moreover, we consider the track direction, and only consider
incidents that are "moving backwards".
"""
# If trains are not travelling in the same direction
# they cannot be due to the same "backward moving" interference
# source.
# (Note: This was the hypothesis when this code was written.
# It turned out that the rogue train could affect all
# trains in the vicinity, not just in the opposite track)
if i["Bound"] != j["Bound"] or \
i["Bound"] not in ['IT', 'OT']:
return False
# time difference in minutes
time_difference = (i["DateTime"] - j["DateTime"]) / np.timedelta64(1, 'm')
location_difference = i["LocID"] - j["LocID"]
if location_difference == 0:
return False
ratio = time_difference / location_difference
if i["Bound"] == 'OT':
return ratio > 0 and ratio < 4
elif i["Bound"] == 'IT':
return ratio < 0 and ratio > -4 [/mw_shl_code]
片段 7
我们发现了满足这个条件的所有事件对：
[mw_shl_code=python,true]incidents = sel_dfincidents.to_records()
# (a, b, c, d, ...) --> ((a,b), (a,c), (a,d), ..., (b,c), (b,d), ..., (c,d), ...)
incident_pairs = list(it.combinations(incidents, 2))
related_pairs = [ip for ip in incident_pairs if same_cascade(*ip)]
related_pairs = [(i[0], j[0]) for i,j in related_pairs] [/mw_shl_code]
片段 8
然后，我们使用并查集的数据结构将所有相关的事件对分组成更大的集合。这使我们可以将可能关联到同一“反常列车”的事件分组。
[mw_shl_code=python,true]def pairs_to_clusters(pairs):
"""
A quick-and-dirty disjoint-set data structure. But this works fast enough for 200+ records
Could be better.
Example input:
(1,2), (2,3), (4,5)
Output:
1: {1,2,3}
2: {1,2,3}
3: {1,2,3}
4: {4,5}
5: {4,5}
"""
the_clusters = dict()
for i,j in pairs:
if i not in the_clusters:
if j in the_clusters:
the_clusters[j].add(i)
the_clusters = the_clusters[j]
else:
the_clusters = set(list([i, j]))
the_clusters[j] = the_clusters
else:
if j in the_clusters:
if the_clusters is not the_clusters[j]: # union the two sets
for k in the_clusters[j]:
the_clusters.add(k)
the_clusters[k] = the_clusters
else: # they are already in the same set
pass
else:
the_clusters.add(j)
the_clusters[j] = the_clusters
return the_clusters [/mw_shl_code]
片段9
然后将我们的算法应用于数据：
[mw_shl_code=python,true]clusters = pairs_to_clusters(related_pairs)
# Show each set only once
clusters = [v for k,v in clusters.items() if min(v) == k]
clusters[0:10] [/mw_shl_code]
片段10
这些是我们确定的一些集群：
[{0, 1},
{2, 4},
{5, 6, 7},
{8, 9},
{18, 19, 20},
{21, 22, 24, 26, 27},
{28, 29, 30, 31, 32, 33, 34},
{42, 44, 45},
{47, 48},
{51, 52, 53, 56}]
接下来，我们计算了可以通过我们的聚类算法解释的事件的百分比。
[mw_shl_code=python,true]# count % of incidents occurring in a cluster
all_clustered_incidents = set()
for i,clust in enumerate(clusters):
all_clustered_incidents |= clust
(len(all_clustered_incidents),
len(incidents),
float(len(all_clustered_incidents)) / len(incidents)) [/mw_shl_code]
片段11
结果是：
[mw_shl_code=python,true](189, 259, 0.7297297297297297) [/mw_shl_code]
它表达的意思是：在数据集中的259个紧急制动事件中，189个案例(其中73个)可以通过“反常列车”假说来解释。我们觉得我们真的走对了路。
我们基于聚类结果对事件图进行了着色。具有相同颜色的三角形在同一个集群中。
图10 通过我们的算法聚类的事件
有多少反常列车?
如图5所示，在环线上的每个端到端行程需要大约1小时。我们通过事件图和与图5密切匹配的线绘制最佳线。这强烈暗示只有一个“反常列车”。
图11 事件集群的时间强烈暗示干扰能够关联到单趟列车
我们还观察到，不明的“反常列车”本身似乎没有遇到任何信号问题，因为它没有出现在我们的散点图。
我们相信我们有一个很好的例子，决定进一步调查。
捕获反常列车
日落之后，我们去了Kim Chuan站以确定“反常列车”。我们不能检查详细的列车日志，因为SMRT需要更多的时间来提取数据。因此，我们决定用传统方式通过查看在事件发生时到达和离开每个车站的列车的视频记录来识别列车。
上午3点，车队出现了头号嫌犯：PV46，一辆自2015年起投入使用的列车。
检验假说
11月6日(星期日)，LTA和SMRT测试，如果PV46是通过运行列车在非高峰时间问题的来源，那我们就对了——PV46确实导致了附近列车之间的通信丢失，并触发那些列车上的紧急制动器。在PV46未投入使用的那天之前，并没有这样的事件发生。
在11月7日(星期一)，我们的团队处理了PV46的历史位置数据，并得出结论，从8月到11月的所有事件中，超过95%可以用我们的假说来解释。其余的事件，可能是由于偶发在正常情况下的信号丢失。
在某些日子，如9月1日，这种模式特别明显。你可以很容易地看到，当PV46投入运行时 ，在其运行期间或者前后，常常发生发生干扰事件。
LTA和SMRT最终于11月11日发布了一份联合新闻稿，与公众分享调查结果。
小结
当我们开始时，我们曾经希望能够找到跨机构调查组可能感兴趣的模式，其中包括LTA、SMRT和DSTA的许多官员。SMRT和LTA提供的整洁的事件日志帮助我们开了一个好头，因为在导入和分析数据之前需要尽可能少的清理。我们还对LTA和DSTA的有效跟踪调查表示满意，因为他们证实了PV46存在硬件问题。
从数据科学的角度来看，我们很幸运，事件发生得如此接近。这使我们能够在如此短的时间内识别问题和罪魁祸首。如果事件更加孤立，则Z字形模式将不那么明显，并且它将让我们耗费更多的时间和数据来解决这个谜。
当然，令人最高兴的是，现在所有人可以坐环线地铁再次充满信心地工作。
注意：本文提到的代码是在2016年11月5日编写的——我们处理SMRT数据以确定环线地铁事件原因的实际日期。我们承认可能存在不足之处。您可以在这里下载我们的Jupyter Notebook的副本(源码)。
来源：大数据杂谈
作者：刘志勇
/2 
京ICP备2020039040号
简书网举报电话：021-34700000
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
千斤顶
显身卡
