熱線電話:13121318867

登錄
首頁精彩閱讀R語言實現交通行業事故案例之黑點確定
R語言實現交通行業事故案例之黑點確定
2016-10-12
收藏

R語言實現交通行業事故案例之黑點確定

淺談道路黑點定義,定義黑點道路為歷史發生事故起數較多和近期發生事故明顯增多兩種道路,并且用簡易事故、一般事故、較大事故、特大事故確定當前發生事故的嚴重程度,即用當量事故數表示,事故越嚴重,則當事事故數越大,當量事故數定義:

一、容易發生黑點道路的確定

1、歷史事故較多道路

通過對各個道路歷史數據的分析,找出歷史發生事故頻率較大的道路作為黑點道路,對于經常發生事故的道路屬于此類。如,取所有道路三年內的當量事故數作為歷史數據,找出當量事故數較大的道路作為預定黑點道路;

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

分析出近期時段較以往事故發生明顯增多道路作為預定黑點道路,這樣可以找出歷史發生事故很少,但是最近明顯發生了很多事故的道路。如,平時最多發生事故起數為1起的事故,近一個月連續發生了3起,則同比增長了200%,則此類道路可作為預定黑點道路。

3、預定黑點道路去重

對1和2分析出的預定黑點道路進行合并,找出所有預定事故黑點道路,因為歷史發生事故較多道路也可能近期突然發生事故數增多,也屬于近期發生事故遽增道路。

二、黑點道路上事故發生較密區域查找

針對確定的預定黑點道路,分別運用聚類算法,找出當前道路上事故發生較密集的各個區域(比如,使用密度聚類算法),作為事故黑點區域。地圖展現時只針對發生較密指定半徑區域為一個事故黑點區(一條道路有可能有個黑點區域),避免地圖展現時整體道路作為一個黑點。

三、黑點位置的確認

根據步驟二分析的事故黑點區域,給定區域中心坐標和半徑在地圖上展現,然后用戶可以標注當前黑點區域的具體位置。

四、具體代碼實現步驟如下

1、連接Oracle數據庫,并讀取所需字段

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

2、分析歷史事故發生較多道路,得到結果集Res

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

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

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

4、預定黑點道路去重,得到結果集Res,并入庫

  1. #3 Combine Res and Res2
  2. DiffRes=Res2[Res2[,1] %in% Res[,1]=='FALSE',]
  3. #B-A
  4. #setdiff(Res2[,1],Res[,1])
  5. Res=rbind(Res,DiffRes)
  6. dbRemoveTable(conn, 'DW_ACD_HDZT1')
  7. 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目的是找到密度相連對象的最大集合。

詳細算法描述參考度娘

  1. #4 DBSCAN(Density-Based Spatial Clustering of Application With Noise)
  2. rs3=dbSendQuery(conn,"select t.lh,t.sgbh,t.sgddzb_x,t.sgddzb_y from ex_dw_acd_file t
  3. join DW_ACD_HDZT1 d on t.lh=d.lh
  4.                 order by lh")
  5. da3=fetch(rs3)
  6. UniqDl3=unique(da3$LH)
  7. da3=data.frame(da3,RCoordinate=as.numeric(da3[,3]),CCoordinate=as.numeric(da3[,4]))
  8. dbRemoveTable(conn, 'DW_ACD_HDZT2')
  9. for(i in 1:length(UniqDl3))
  10. {
  11.  #i=1
  12.  tempdate=0
  13.  tempdate=da3[da3[,1]==UniqDl3[i],]
  14.  #Normalization
  15.  tempdate[,5]=(tempdate[,5]-min(tempdate[,5]))/(max(tempdate[,5])-min(tempdate[,5]))
  16.  tempdate[,6]=(tempdate[,6]-min(tempdate[,6]))/(max(tempdate[,6])-min(tempdate[,6]))
  17.  #version one
  18.  #set.seed(664554)
  19.  #n <- 10
  20.  #x <- cbind(runif(10, 0, 10)+rnorm(n, sd=0.4), runif(10, 0, 10)+rnorm(n,sd=0.4))

  21.  #help(dbscan)
  22.  #ds <- dbscan(x, 0.05,MinPts=5,showplot=1)
  23.  ds <- dbscan(tempdate[,5:6], 0.005,MinPts=5,showplot=1) # run with showplot=1 to see how dbscan works.
  24.  aa=as.matrix(ds$cluster)
  25.  RE=as.data.frame(cbind(tempdate,aa))
  26.  colnames(RE)=c( "LH","SGBH","SGDDZB_X","SGDDZB_Y","RCoordinate","CCoordinate","ClusterNum")
  27.  dbWriteTable(conn,'DW_ACD_HDZT2',RE, row.names = F, append = TRUE)
  28. }
  29. #seq.POSIXt(ISOdate(2001,1,1)-60*60*12,ISOdate(2001,1,3)+60*60*12,"hours")


數據分析咨詢請掃描二維碼

若不方便掃碼,搜微信號:CDAshujufenxi

數據分析師資訊
更多

OK
客服在線
立即咨詢
日韩人妻系列无码专区视频,先锋高清无码,无码免费视欧非,国精产品一区一区三区无码
客服在線
立即咨詢