名侦探R语言，巧破交通事故黑点

原创 2016-11-23 华青莲 R语言中文社区

简述：伴随着人们对出行以及货物运输的安全、快速、方便和机动性要求的提高，公路交通事故及其造成的损失越来越受世界各国高度重视,道路交通安全已成为社会各界广泛关注的重大课题。为了减少公路交通事故发生次数、降低事故严重程度，人们一直在进行着不懈的努力和探索。公路交通事故黑点分析是公路交通安全研究中重要的内容之一，鉴别出公路交通事故黑点并分析其存在的主要原因，从而提出切实可行的交通治理措施以改善公路上的交通安全状况，是公路交通管理工作中的重要工作内容，也是社会价值观、科学研究应用、政策制定实施中”以人为本“思想的直接体现。浅谈道路黑点定义，定义黑点道路为历史发生事故起数较多和近期发生事故明显增多两种道路，并且用简易事故、一般事故、较大事故、特大事故确定当前发生事故的严重程度，即用当量事故数表示，事故越严重，则当事事故数越大，当量事故数定义：

一、容易发生黑点道路的确定

1、历史事故较多道路

通过对各个道路历史数据的分析，找出历史发生事故频率较大的道路作为黑点道路，对于经常发生事故的道路属于此类。如，取所有道路三年内的当量事故数作为历史数据，找出当量事故数较大的道路作为预定黑点道路；

2、近期发生事故遽增道路

分析出近期时段较以往事故发生明显增多道路作为预定黑点道路，这样可以找出历史发生事故很少，但是最近明显发生了很多事故的道路。如，平时最多发生事故起数为1起的事故，近一个月连续发生了3起，则同比增长了200%，则此类道路可作为预定黑点道路。

3、预定黑点道路去重

对1和2分析出的预定黑点道路进行合并，找出所有预定事故黑点道路，因为历史发生事故较多道路也可能近期突然发生事故数增多，也属于近期发生事故遽增道路。

二、黑点道路上事故发生较密区域查找

针对确定的预定黑点道路，分别运用聚类算法，找出当前道路上事故发生较密集的各个区域（比如，使用密度聚类算法），作为事故黑点区域。地图展现时只针对发生较密指定半径区域为一个事故黑点区（一条道路有可能有个黑点区域），避免地图展现时整体道路作为一个黑点。

三、黑点位置的确认

根据步骤二分析的事故黑点区域，给定区域中心坐标和半径在地图上展现，然后用户可以标注当前黑点区域的具体位置。

具体代码实现步骤如下：

1、连接Oracle数据库，并读取所需字段

library(DBI)
library(ROracle)
library(cluster)
#install.packages("fpc")
library(fpc)
#w=c(1,0.5,0.3,0.1,0.01)
drv=dbDriver('Oracle')
conn=dbConnect(drv,'acci_tz','acci_tz','localhost:1521/jgyw')

2、分析历史事故发生较多道路，得到结果集Res

#1 Get the history black spots
rs=dbSendQuery(conn,"select lh,sum(dlsgs) dlsgs from ex_dw_acd_file
 group by lh order by lh asc")
da=fetch(rs)
#ps=data.frame(paste(da[,1],da[,2],sep=''),da[,3])  #splice
#test
class<-c(1,2,3,2,1,2,1,3)
c(81,65,72,88,73,91,56,90)->ca
factor(class)->class
tapply(ca, class, mean)
# divide into groups
ps1=as.data.frame(tapply(da[,2],da[,1],sum))    
#Get the frequency of each road and cumulative frequency
ps2=as.data.frame(ftable(ps1[,1]))
ps3=cumsum(ps2[,2]/(sum(ps2[,2])))
ps3=data.frame(ps2[,1],ps3)
#plot(ps3)
#points(ps3,type='l')
 #Set the Threshold
yz=ps3[ps3[,2]>=0.95,1][1] 
 #Black Spots Judgement
ps4=data.frame(ps1[ps1[,1]>=as.numeric(as.character(yz)),])
# Get In 
ZT=as.data.frame(as.integer(da[,1] %in%  row.names(ps4))) 
c=data.frame(da[,1],ZT)
Res=data.frame(c[c[,2]==1,])
colnames(Res)<-c('LH','ZT')

3、分析近期发生事故遽增道路Res2

#2 Get the recent black spots
drv=dbDriver('Oracle')
conn=dbConnect(drv,'acci_tz','acci_tz','localhost:1521/jgyw')
rs2=dbSendQuery(conn,"select substr(d.sgfssj,1,6) yfbm,d.lh,sum(d.dlsgs) dlsgs 
from ex_dw_acd_file d where 
                substr(d.sgfssj,1,6) between 201511 and 201601
                group by substr(d.sgfssj,1,6),d.lh
                order by d.lh,substr(d.sgfssj,1,6) asc")
da2=fetch(rs2)
# To be unique of road
UniqDl=unique(da2$LH)
Res2=c()
for (i in 1:length(UniqDl))
{
  #i=15
  tempd=da2[da2[,2] %in% UniqDl[i],]
  #Filter the three months data
  if(nrow(tempd)>2)
  {
    #Filter the  Slope of difference which is more then 1.5
    if(tempd[3,3]/tempd[1,3]>=1.5 && tempd[3,3]/tempd[2,3]>=1.5)
    {
      Res2=rbind(Res2,UniqDl[i])
    }
  }
}
Res2=cbind(Res2,rep(1,length(Res2)))
Res2=as.data.frame(Res2)
colnames(Res2)<-c('LH','ZT')

4、预定黑点道路去重，得到结果集Res，并入库

#3 Combine Res and Res2
DiffRes=Res2[Res2[,1] %in% Res[,1]=='FALSE',]
#B-A
#setdiff(Res2[,1],Res[,1])
Res=rbind(Res,DiffRes)
dbRemoveTable(conn, 'DW_ACD_HDZT1')
dbWriteTable(conn,'DW_ACD_HDZT1',Res, row.names = F, append = TRUE)

5、黑点道路上事故发生较密区域查找，使用密度聚类算法DBSCAN

附DBSCAN：

DBSCAN(Density-Based Spatial Clustering of Applications with Noise)是一个比较有代表性的基于密度的聚类算法。与划分和层次聚类方法不同，它将簇定义为密度相连的点的最大集合，能够把具有足够高密度的区域划分为簇，并可在噪声的空间数据库中发现任意形状的聚类。DBSCAN自动地确定簇个数，而对于K-means，簇个数需要作为参数指定。然而，DBSCAN必须指定另外两个参数：Eps（邻域半径）和MinPts（最少点数）。

DBSCAN中的几个定义：

Ε邻域：

给定对象半径为Ε内的区域称为该对象的Ε邻域；

核心对象：如果给定对象Ε领域内的样本点数大于等于MinPts，则称该对象为核心对象；

直接密度可达：

对于样本集合D，如果样本点q在p的Ε领域内，并且p为核心对象，那么对象q从对象p直接密度可达。

密度可达：

对于样本集合D，给定一串样本点p1,p2….pn，p= p1,q= pn,假如对象pi从pi-1直接密度可达，那么对象q从对象p密度可达。

密度相连：

存在样本集合D中的一点o，如果对象o到对象p和对象q都是密度可达的，那么p和q密度相联。

可以发现，密度可达是直接密度可达的传递闭包，并且这种关系是非对称的。密度相连是对称关系。DBSCAN目的是找到密度相连对象的最大集合。

详细算法描述参考度娘

#4 DBSCAN(Density-Based Spatial Clustering of Application With Noise)
rs3=dbSendQuery(conn,"select t.lh,t.sgbh,t.sgddzb_x,t.sgddzb_y from ex_dw_acd_file t
join DW_ACD_HDZT1 d on t.lh=d.lh
                order by lh")
da3=fetch(rs3)
UniqDl3=unique(da3$LH)
da3=data.frame(da3,RCoordinate=as.numeric(da3[,3]),CCoordinate=as.numeric(da3[,4]))
dbRemoveTable(conn, 'DW_ACD_HDZT2')
for(i in 1:length(UniqDl3))
{
 #i=1
 tempdate=0
 tempdate=da3[da3[,1]==UniqDl3[i],] 
 #Normalization
 tempdate[,5]=(tempdate[,5]-min(tempdate[,5]))/(max(tempdate[,5])-min(tempdate[,5]))
 tempdate[,6]=(tempdate[,6]-min(tempdate[,6]))/(max(tempdate[,6])-min(tempdate[,6]))
 #version one
 #set.seed(664554)
 #n <- 10
 #x <- cbind(runif(10, 0, 10)+rnorm(n, sd=0.4), runif(10, 0, 10)+rnorm(n,sd=0.4))
 
 #help(dbscan)
 #ds <- dbscan(x, 0.05,MinPts=5,showplot=1)
 ds <- dbscan(tempdate[,5:6], 0.005,MinPts=5,showplot=1) # run with showplot=1 to see how dbscan works.
 aa=as.matrix(ds$cluster)
 RE=as.data.frame(cbind(tempdate,aa)) 
 colnames(RE)=c( "LH","SGBH","SGDDZB_X","SGDDZB_Y","RCoordinate","CCoordinate","ClusterNum")
 dbWriteTable(conn,'DW_ACD_HDZT2',RE, row.names = F, append = TRUE)
}
#seq.POSIXt(ISOdate(2001,1,1)-60*60*12,ISOdate(2001,1,3)+60*60*12,"hours")