한국 지진 데이터의 군집 모델

개요

데이터는 고감자 (gogamza)님의 누리집 http://freesearch.pe.kr/archives/2884에서 참고 했다. 이 사이트는 ggmap과 지진 강도 추이와 이상점(outlier)에 대한 시각화가 잘 되어 있다. 여기서는 각종 데이터 군집 모델을 수행하고 결과를 다른 정보(한국 지질 구조)와 비교한다.

1. 데이터 준비

고감자님의 데이터 주소에서 다운 받아 csv로 저장한 후 데이터 전처리 작업을 수행한다. 먼저 데이터를 확인한다.

> data.tmp <- read.csv2(file="./data/eq.csv", header=T, stringsAsFactors = F, sep = "\t")
> head(data.tmp)
              date power latitude longitude                               place
1  2012-06-24 1:20   2.1  37.53 N  124.56 E      인천 백령도 남남서쪽 48km 해역
2 2012-06-23 16:42   2.6   36.7 N  128.31 E        경북 문경시 북동쪽 17km 지역
3 2012-06-19 20:55   3.1  33.48 N  125.93 E   제주 제주시 고산 북서쪽 30km 해역
4  2012-06-14 6:15   2.1  39.68 N  126.81 E        함경남도 요덕 북서쪽 5km지역
5 2012-06-12 17:52   2.2  35.87 N  128.67 E       대구 수성구 동북동쪽 4km 지역
6 2012-06-12 13:34   2.4  37.33 N  125.27 E 인천 옹진군 연평도 남서쪽 53km 해역

데이터 정제 작업을 수행한다.

# 위도와 경도를 수치로 표현하기 위한 함수
> mapCoordinate2number <- function (str){
    flag <- substr(str, nchar(str), nchar(str))
    num <- substr(str, 1, nchar(str) - 1)
    
    num <- as.numeric(num)
    if (flag =="W"  || flag =="S" )
      num = -num
      
    return (num)
}
  
> data.eq <- transform(data.tmp, 
                     date = as.POSIXlt.date(strftime(date, format = "")),
                     latitude = mapCoordinate2number(latitude),
                     longitude = mapCoordinate2number(longitude),
                     power = as.numeric(power))

> str(data.eq)
'data.frame': 978 obs. of  5 variables:
 $ date     : POSIXlt, format: "2012-06-24 01:20:00" "2012-06-23 16:42:00" "2012-06-19 20:55:00"  ...
 $ power    : num  2.1 2.6 3.1 2.1 2.2 2.4 2.1 2.1 2.1 2.6 ...
 $ latitude : num  37.5 36.7 33.5 39.7 35.9 ...
 $ longitude: num  125 128 126 127 129 ...
 $ place    : chr  "인천 백령도 남남서쪽 48km 해역" "경북 문경시 북동쪽 17km 지역"  ...

데이터 요약 정보를 혹인하고 결측값(NA)가 있는 행(row)를 제거한다.

> summary(data.eq[2:4])
     power          latitude       longitude    
 Min.   :1.700   Min.   :32.35   Min.   :122.8  
 1st Qu.:2.400   1st Qu.:35.76   1st Qu.:125.9  
 Median :2.700   Median :36.58   Median :127.1  
 Mean   :2.769   Mean   :36.66   Mean   :127.2  
 3rd Qu.:3.000   3rd Qu.:37.80   3rd Qu.:128.5  
 Max.   :5.300   Max.   :41.60   Max.   :131.1  
                 NA's   :2       NA's   :2      

> data.eq <- data.eq[complete.cases(data.eq[,3:4]),]
> summary(data.eq[2:4])
     power          latitude       longitude    
 Min.   :1.700   Min.   :32.35   Min.   :122.8  
 1st Qu.:2.400   1st Qu.:35.76   1st Qu.:125.9  
 Median :2.700   Median :36.58   Median :127.1  
 Mean   :2.769   Mean   :36.66   Mean   :127.2  
 3rd Qu.:3.000   3rd Qu.:37.80   3rd Qu.:128.5  
 Max.   :5.300   Max.   :41.60   Max.   :131.1 

구글 지도(ggmap)를 이용하여 지진 데이터의 분포를 확인한다. 육안으로 결과 지도를 확인해도 별다른 특징을 잡아내기 힘들다는 것을 알게 될 것이다. 단지 몇 개의 군집을 이루고 있다는 것을 알 수 있다.

> library(ggmap)

> map.hdf1 <- get_map("seoul, korea", zoom = 5, maptype = "terrain")  #대한민국
> map.hdf2 <- get_map("seoul, korea", zoom = 6, maptype = "terrain")  #남한
> map.hdf3 <- get_map("seoul, korea", zoom = 7, maptype = "terrain")  #남한

# Figure 1를 그림을 출력
> ggmap(map.hdf2) + 
    geom_point(aes(x = longitude, y = latitude, colour =power), data = data.eq, alpha = .5) + 
    scale_colour_gradient(low = "lightblue", high = "blue")

Figure 1: 한국 지진 데이터의 시작화

구글맵을 이용하여 한국 지진 데이터를 시각화 했지만 바로 눈에 들어오는 특징이 없다. 이 데이터를 군집화 하면 더 심도 있게 이해할 수 있다. 다음 장부터는 군집 모델을 하나씩 사용한다. 일반적으로 데이터 군집은 데이터 특징을 이해하고 다른 학습 모델을 적용하기 위한 발판으로 이용된다.

2. 군집(Clustering)

2.1. The k-Means Clustering

이 절에서는 k-means 군집 모델을 수행한다. 군집은 할당된 중앙점과 가장 가까운 제곱합과 같은 k개의 그룹으로 분할하는 것이 목적이다. 모델을 더 이해하기 위해 algorithm와 method 변수를 조정하여 모델을 실행해야 한다.

우리는 지진 데이터의 특징을 알 지 못하기 때문에 군집 개수를 변경하면서 육안으로 가장 적절하게 보이는 군집의 개수를 선택하여 다른 군집 모델과 비교한다.

이전 절에서 사용한 ggmap의 결과를 그래프를 계속 사용한다.

# -------------------------------
# 20개의 군집 생성
kmenas.20 <- kmeans(data.eq[3:4], 20)

# 이전 장에 사용한 ggmap을 사용
map.eq <- ggmap(map.hdf2) + geom_point(aes(x = longitude, y = latitude, color= power), data = data.eq) +
  scale_colour_gradient(low = "lightblue", high = "blue")


# 군집을 결합
map.eq.cluser <- map.eq + geom_point(aes(x = longitude, y = latitude), 
                                     data = as.data.frame(kmenas.20$centers), 
                                     color="red", fill = "red",  size = 2, show.legend = F)


#화면 출력
print(map.eq.cluser)


# -----------------------------------------------
# 40개의 군집 생성
kmenas.40 <- kmeans(data.eq[3:4], 40)
map.eq.cluser <- map.eq + geom_point(aes(x = longitude, y = latitude), 
                                     data = as.data.frame(kmenas.40$centers), 
                                     color="black", fill = "black",  size = 2, show.legend = F)

print(map.eq.cluser)

# -----------------------------------------------
# 60개의 군집 생성
kmenas.60 <- kmeans(data.eq[3:4], 60)
map.eq.cluser <- map.eq + geom_point(aes(x = longitude, y = latitude), 
                                     data = as.data.frame(kmenas.60$centers), 
                                     color="black", fill = "black",  size = 2, show.legend = F)

print(map.eq.cluser)


# -----------------------------------------------
# 80개의 군집 생성
kmenas.80 <- kmeans(data.eq[3:4], 80)

map.eq.cluser <- map.eq + geom_point(aes(x = longitude, y = latitude), 
                                     data = as.data.frame(kmenas.80$centers), 
                                     color="black", fill = "black",  size = 2, show.legend = F)

Figure 2: k-means모델에서 20개의 군집

Figure 3: k-means모델에서 40개의 군집

Figure 4: k-means모델에서 60개의 군집

Figure 5: k-means모델에서 80개의 군집

결과 그림을 클릭해서 확대해 보면 군집의 수가 증가할수록 어떤 선이 보이는 경향이 강하다. 이 선이 어떤 선을 표현하는지는 다른 정보와 비교하면 더 명확하게 이해할 수 있다.

2.2. The k-Medoids Clustering

이 절은 pam()과 pamk()로 수행하는 k-medoids 모델을 보여준다. K-means가 군집의 중앙점과 가장 가까운 곳에 할당되지만, k-medoids는 군집의 중앙과 가장 가까운 객체를 가지고 표현한다. 이상점에 더 강건한 특징을 가진다. 여기서는 군집의 개수를 결정할 수 있는 pamk()를 사용하여 모델링한다.

# pam, pamk를 수행하는 패키지
library(fpc)

# 40개의 군집 생성
pam.40 <- fpc::pamk(data.eq[,3:4], 40)

map.eq.pamk.40 <- map.eq + geom_point(aes(x = longitude, y = latitude), 
                                     data = as.data.frame(pam.40$pamobject$medoids), 
                                     color="black", fill = "black",  size = 2, 
                                     show.legend = F)

print(map.eq.pamk.40)

Figure 3: k-means모델에서 40개의 군집

Figure 6: k-medoids모델에서 40개의 군집

그림에서 k-means 와 k-medoids 의 결과는 차이가 거의 없음을 알 수 있다. 단지 중앙점의 위치가 약간 다르다.

2.3. The k-Medoids Clustering

이 함수는 n개의 객체에서 비유사도(dissimilarities)의 집합을 사용하여 hierarchical clustering 분석을 수행한다.
여기서는 약간의 코딩이 필요하다

hc <- hclust(dist(data.eq[,3:4]), method="ave")

# 트리 구조의 그래프를 그림, 그림 없음
plot(hc, hang = -1, labels=F)
rect.hclust(hc, k=40)
groups <- cutree(hc, k=40)

# 군집의 평균 구함
hc.40 <- aggregate(data.eq[, 3:4], by = list(group = tmp$group), FUN = mean)

# 지도 그래프 그림
map.eq.hc.40 <- map.eq + geom_point(aes(x = longitude, y = latitude), 
                                      data = hc.40, 
                                      color="black", fill = "black",  size = 2,
                                      show.legend = F)

print(map.eq.hc.40)

Figure 7: hierarchical clustering 모델에서 40개의 군집

2.4. Density-based Clustering

밀도 기반 군집의 사상은 만약 밀집 분포 지역의 다른 하나와 연결된다면 군집 객체를 같은 클러스터에 속해야 한다.
여기에 중요한 변수 2개가 있다

• eps : 이웃의 크기를 결정하는 도달 거리
• MinPts : 최소 군집 수

library(fpc)

# eps, MinPts 의 값을 조절하면서 모델링 해야 함
clstr.ds <- dbscan(data.eq[,3:4], method = "hybrid", eps=0.2, MinPts=5)
ds.40 <- aggregate(data.eq[, 3:4], by = list(group = clstr.ds$cluster), FUN = mean)

map.eq.ds <- ggmap(map.hdf2) + 
    geom_point(aes(x = longitude, y = latitude, color= clstr.ds$cluster, 
                   legend="cluster"), 
               data = data.eq) + 
    scale_colour_gradient(low = "lightgreen", high = "darkgreen")
  

map.eq.ds.40 <- map.eq.ds + geom_point(aes(x = longitude, y = latitude), 
                                    data = ds.40 , 
                                    color="black", fill = "black",  size = 2,
                                    show.legend = F)

map.eq.ds.40 <- map.eq.ds.40 +labs(colour="cluster")


print(map.eq.ds.40)

Figure 8: dbscan() 함수에의한 밀도 기반 군집

2.5. Kohonet 패키지의 자가 조직화 맵(som)

som 모델은 Euclidean 거리 계산 공식을 사용하여 구현되었다. som 구현에 있어 군집의 개수를 결정은 somgrid 함수에서 중요 변수 xdim, ydim 값을 설정함으로써 이루어진다. 기본값 xdim = 8, ydim = 6 이다. 기본 설정으로 som 모델을 구현하면 48개의 군집을 이룬다.

library(kohonen)

data.map <- data.eq[, c("latitude", "longitude")]

# 군집의 개수를 결정하기 위한 차원 설정, 8 X 6 = 48 개의 군집이 생성됨
som_grid <- somgrid(xdim = 8, ydim=6, topo="hexagonal")

som_model <- kohonen::som(as.matrix(data.map), 
                 grid=som_grid, 
                 keep.data = TRUE,
                 n.hood="circular" )

# summary(som_model)
# print(som_model)

# 지도에서 사용할 데이터 생성
data.map.somgroups <- cbind(data.map, classif = som_model$unit.classif)

# 군집을 표현할 데이터
tmp <- aggregate(data.map, by = list(classif = som_model$unit.classif), mean)

library(ggmap)
library(ggplot2)

hdf2 <- get_map("seoul, korea", zoom = 6, maptype = "terrain")  #남한

map <- ggmap(hdf2) + 
  geom_point(aes(x = longitude, y = latitude,  group=classif, color=factor(classif), 
                 fill=factor(classif)), 
             data = data.map.somgroups, alpha = .5) +
  geom_point(data = tmp, aes(x = longitude, y = latitude,  shape= "."), 
             show.legend = F, size=3)

# Figure 9 그림 
plot(map)

Figure 9: som 모델로 생성한 군집

2.6. randomUnfiformForest를 이용한 군집

randomUniformForest를 이용하여 클러스터 작업을 수행한다.

library(randomUniformForest)

data.map <- data.eq[, c("latitude", "longitude")]

# endModel 기본, 클러스터를 50개로 설정
fit.ruf = unsupervised.randomUniformForest(data.map, 
                                           endModel = "MDSkMeans", # 기본값
                                           clusters = 50)

# 지도에 표현할 데이터 생성
data.map.ruf.cluster <- cbind(data.map, classif = fit.ruf$unsupervisedModel$cluster)

# 클러스터 중심 좌표 생성
tmp <- aggregate(data.map, by = list(classif = fit.ruf$unsupervisedModel$cluster), mean)

library(ggmap)
library(ggplot2)

hdf2 <- get_map("seoul, korea", zoom = 6, maptype = "terrain")  #남한

map <- ggmap(hdf2) + 
  geom_point(aes(x = longitude, y = latitude,  group=classif, color=factor(classif), 
                 fill=factor(classif)), 
             data = data.map.ruf.cluster, alpha = .5) +
  geom_point(data = tmp, aes(x = longitude, y = latitude,  shape= "."), 
             show.legend = F, size=2)

# Figure 10
plot(map)

Figure 10: endModel = "MDSkMeans" 기본값으로 클러스터

randomUniformForest의 endModel 변수를 SpectralkMeans로 설정하여 클러스터 작업을 수행한다.

library(randomUniformForest)

data.map <- data.eq[, c("latitude", "longitude")]

fit.ruf.spectral = unsupervised.randomUniformForest(data.map, 
                                                    endModel = "SpectralkMeans",
                                                    clusters = 50)


data.map.ruf.cluster <- cbind(data.map, 
                              classif = fit.ruf.spectral$unsupervisedModel$cluster)
tmp <- aggregate(data.map, 
                 by = list(classif = fit.ruf.spectral$unsupervisedModel$cluster), mean)

library(ggmap)
library(ggplot2)

hdf2 <- get_map("seoul, korea", zoom = 6, maptype = "terrain")  #남한

map <- ggmap(hdf2) + 
  geom_point(aes(x = longitude, y = latitude,  group=classif, color=factor(classif), 
                 fill=factor(classif)), 
             data = data.map.ruf.cluster, alpha = .5) +
  geom_point(data = tmp, aes(x = longitude, y = latitude,  shape= "."), 
             show.legend = F, size=2)

# Figure 11
plot(map)

Figure 11: endModel = "SpectralkMeans" 기본값으로 클러스터. 클러스터 중심이 뚜렷하게 한반도를 둘러싸고 있음

3. 한반도 지각 구조와 군집 데이터의 비교

한반도의 지각 구조와 군집결과를 비교해 보면 어느 정도 관련성이 있음이 나타난다. K-means에서 80개의 군집그립과 비교해 보라. 군집의 개수가 증가할수록 지각의 분할되는 선과 더 비슷하게 군집의 중앙점들이 나타나는 것을 확인할 수 있다. 다만 너무 군집의 수가 증가하면 더 불일치하게 보인다.

Figure 12: 한반도 지각 구조

한반도는 유라시아판에 동남쪽 가장자리에 있다. 이쪽에 영향을 주는 힘은 동쪽에서 태평양 판 그리고 남쪽에서 필리핀 판이다. 이 힘을 한반도 지각에 영향을 받는다고 가정하면 지진의 발생 위치가 태백산맥 줄기를 기준으로 동, 남, 서쪽 가장자리에서 잦은 빈도로 발생하는 점이 이해할 수 있다. 자세히 정리된 블로그 주소는 http://www.seehint.com/r.asp?no=12937 이다.

Figure 13: GPS로 측정한 지구 지각의 운동 방향과 속도