前一篇文章講到了軌跡的繪制，那么今天就順著講下軌跡的聚類，比如說我們常見的，寒潮研究把寒潮的冷空氣路徑分為三類(西路，北路，西北路)，或者臺風(fēng)研究中也可以根據(jù)臺風(fēng)的速度，或者拐點(diǎn)等要素對臺風(fēng)軌跡分類研究，再比如降水可以對不同降水事件的水汽來源做聚類研究等等。
實(shí)際上NCL就可以實(shí)現(xiàn)聚類算法，官方提供了k-means的算法函數(shù)(kmeans_as136)，這個函數(shù)筆者也用過，相比與Python的聚類可以說還是比較麻煩的。
另外我也見過大佬們用MetoInfo的軟件做軌跡和聚類研究，十分方便。
K-means聚類
百度搜索python的K-means算法，可以得到的結(jié)果還是很多的，然而大部分都是基于鳶尾花數(shù)據(jù)實(shí)現(xiàn)的，例子有限，那我就結(jié)合我所使用的寒潮路徑數(shù)據(jù)來做個印證，對比起來更容易理解。
下圖是我們的總路徑，我們現(xiàn)要對其進(jìn)行分類。
使用的數(shù)據(jù)是這樣的：
print(x.shape)
print(y.shape)
#(216, 68)，216為216條路徑，68指每條路徑由68個時刻點(diǎn)構(gòu)成。
#(216, 68)
(這套數(shù)據(jù)是使用FLEXPART后向追蹤得到的，與我目前的論文有關(guān)，因此暫時不提供)
首先調(diào)用庫
from sklearn.cluster import KMeans
接著整合數(shù)據(jù)，因?yàn)槲覀兊腦和Y是兩個數(shù)組存放的，現(xiàn)在合并起來，這一步的方法有很多(concatenate,stack等等都可以)，我給一個容易理解的方法：
traj = np.zeros((216,68,2))
traj[:,:,0] = x
traj[:,:,1] = y
traj = traj.reshape((216,68*2))
最后的reshape我們得到了一個[216,136]的數(shù)組，要知道聚類是不會管你的點(diǎn)的實(shí)際物理意義的，對他來說都是數(shù)字而已，每一個時刻的一個點(diǎn)是由2個要素構(gòu)成的，即經(jīng)度和維度，那么1條完整的路徑是由68*2個要素構(gòu)成的，因此我們?nèi)ゾ垲悤r，只要將需要考慮的要素全部按順序給進(jìn)去就可以了，但是有一個問題是需要注意的，聚類時，如果要素過大，而樣本過少，就會導(dǎo)致信息熵很大，這樣的聚類雖然不能說有錯，但是起碼誤差很大，能否達(dá)到需要的效果就需要慎重考慮了，因此，當(dāng)要素過大時，可以考慮用PCA方法或者其他手段，提取主成分，降低信息熵，對于一般的軌跡聚類，是不用考慮這些的，就比如這個例子，雖然我們給出了136個要素，但是效果還是符合預(yù)期的。
那么接下來就可以聚類了
km = KMeans(n_clusters=3)#構(gòu)造聚類器
km.fit(traj)#聚類
label = km.labels_ #獲取聚類標(biāo)簽
K-means需要我們實(shí)現(xiàn)給定標(biāo)簽，這與層次聚類不同，層次聚類是兩兩最優(yōu)合并，因此會得出一個聚類樹狀圖，根據(jù)結(jié)果選定聚類數(shù)。K-means則是事先給定聚類數(shù)，然后計(jì)算聚類效果(通過一些參數(shù)評估)，不斷調(diào)整聚類數(shù)，最終確定聚類數(shù)。
可以看到，216條路徑已經(jīng)被分為了3種，我們現(xiàn)在逐類挑出即可
label0=np.array(np.where(labels==0))
label1=np.array(np.where(labels==1))
label2=np.array(np.where(labels==2))
numlabel0=len(label0[0,:])
numlabel1=len(label1[0,:])
numlabel2=len(label2[0,:])
#[0,:]是因?yàn)樯厦嫒〕鰜淼臄?shù)據(jù)是兩維的，第一維是1，因此需要指定才可以獲取長度
x0 = np.array([x[i,:] for i in label0]).reshape((numlabel0,68))
y0 = np.array([y[i,:] for i in label0]).reshape((numlabel0,68))
x1 = np.array([x[i,:] for i in label1]).reshape((numlabel1,68))
y1 = np.array([y[i,:] for i in label1]).reshape((numlabel1,68))
x2 = np.array([x[i,:] for i in label2]).reshape((numlabel2,68))
y2 = np.array([y[i,:] for i in label2]).reshape((numlabel2,68))
這段就可以看出，python的代碼還是很簡潔的。
那現(xiàn)在我們獲取了每一類的經(jīng)度和緯度，按我上一篇文章的方法繪制即可。
K-means聚類效果評估
最后，講一下K-means聚類的效果如何評估，通常用幾個參數(shù)來確定聚類數(shù)。這個要用數(shù)據(jù)說話，不能很主觀的我們想分幾類就分幾類。
一、SSE（簇內(nèi)誤方差）
SSE參數(shù)的核心思想就是計(jì)算誤方差和，SSE的值越小，證明聚類效果越好，當(dāng)然，聚類數(shù)越大，SSE必然是越小的，SSE的分布類似對數(shù)函數(shù)，是逐漸趨0的，同樣也類似對數(shù)函數(shù)，有一個突然的拐點(diǎn)，即存在一個下降趨勢突然變緩的點(diǎn)，這個點(diǎn)對應(yīng)的K值即為最佳聚類數(shù)。
SSE = []
for i in range(1,11):
km = KMeans(n_clusters=i)
km.fit(traj)
#獲取K-means算法的SSE
SSE.append(km.inertia_)
可以看到，k=3之后的誤方差和下降變緩，因此3類是最佳選擇。
二、輪廓系數(shù)S
對于其中的一個點(diǎn) i 來說：
計(jì)算 a(i) 為i向量到所有它屬于的簇中其它點(diǎn)的距離的平均值
計(jì)算 b(i) i向量到各個非本身所在簇的所有點(diǎn)的平均距離的最小值
那么 i 向量輪廓系數(shù)就為：S(i) = (b(i)-a(i))/(max(a(i),b(i)))
當(dāng)然了，公式并不重要，如果需要在論文中給出的話，建議查詢維基百科。
SKLEANR庫給出了計(jì)算函數(shù)。
from sklearn.metrics import silhouette_score
S = [] # 存放輪廓系數(shù)
for i in range(2,10):
kmeans = KMeans(n_clusters=i) # 構(gòu)造聚類器
kmeans.fit(traj)
S.append(silhouette_score(traj,kmeans.labels_,metric='euclidean'))
需要注意的是，這里使用的衡量標(biāo)準(zhǔn)是"euclidean"也就是常說的歐氏距離。
可以看到，k=3時，局部輪廓系數(shù)最大(極大值)，因此K=3為最優(yōu)選擇。
如何確定數(shù)據(jù)是否合適使用K-means方法
通常的話前面介紹到的兩種評估參數(shù)即可確定該數(shù)據(jù)能否使用K-means聚類方法，不需要兩種同時使用，只需要其中一種評估參數(shù)可以挑選出一個最優(yōu)K值即可。如果兩種參數(shù)都無法挑出最優(yōu)K值，那么說明你可能需要另外再找一種聚類方法了，再或者就是我前文提到的，數(shù)據(jù)要素太大，信息熵溢出了，使用PCA提取主成分，對數(shù)據(jù)進(jìn)行降維。
