1. 為什么Kmeans對原始軌跡坐標“失效”了三條最容易犯的錯大概兩年前我接過一個共享電單車的騎行軌跡分析需求客戶只說了一句“我們用Kmeans聚一下就行”。我當時心里就咯噔一下因為軌跡聚類這個事最容易翻車的地方恰恰就在“直接聚”這三個字上。后來果不其然他們對原始GPS坐標點跑了Kmeans出來的圖完全看不出任何騎行模式每條軌跡都像被剪刀剪碎了一樣東一塊西一塊。這件事之后我養(yǎng)成了一個習慣拿到軌跡數(shù)據(jù)的第一反應不是寫Kmeans而是先想清楚“我要聚的對象到底是什么”。Kmeans本身只是一個在歐氏空間里找簇心的算法它天生處理的是“點”和“點”之間的距離。但軌跡是一條帶順序和時間的序列同樣是從A點到B點有人走直線有人繞路有人在某個路口等了兩分鐘紅燈有人中途停下去買東西。這些行為差異如果直接反映在坐標點上Kmeans根本區(qū)分不出“這是一條繞路的軌跡”和“這是兩個不同地點的散點”。所以這篇文章要解決的問題就是怎樣把Kmeans正確用在軌跡數(shù)據(jù)上并且給出可以直接跑的Matlab代碼。我把新手最容易犯的錯總結(jié)成三條先看明白這三條后面再談算法路線和代碼思路會順很多。1.1 軌跡長度不一致特征矩陣根本構(gòu)造不起來這是最樸素也最致命的問題。Kmeans的輸入是一個二維矩陣每一行是一個樣本每一列是一個特征。但軌跡數(shù)據(jù)天生長短不一A用戶騎了5分鐘GPS采了30個點B用戶騎了40分鐘采了300個點。你沒法直接把這兩條軌跡放進行矩陣因為維度都對不上。大多數(shù)人的第一反應是“補零”或“截斷到最短”。補零的問題在于零本身會成為一個有意義的特征Kmeans會莫名其妙多分出一類專門吸收那些被補了超多零的短軌跡截斷更粗暴直接把軌跡后半段的形態(tài)扔掉了比如那條繞路去充電樁換電的軌跡后半段恰恰是它和其他軌跡最大的區(qū)別。正確做法是先對軌跡做弧長重采樣把所有軌跡統(tǒng)一成固定數(shù)量的坐標點。這一步做好了后續(xù)的Kmeans才有“矩陣”可吃。1.2 把所有坐標點堆在一起聚類出來的不是軌跡簇而是熱力點這是“假聚類”里面最多人踩的坑。有些人知道軌跡長度不一但嫌重采樣麻煩索性把100條軌跡的上萬個點全部倒進同一個矩陣然后跑Kmeans。跑完之后聚類結(jié)果確實很干凈簇與簇之間邊界明顯但每個簇的含義是“這些坐標點在空間上挨得比較近”不是“這些軌跡在行為模式上相似”。舉個很簡單的例子三條軌跡都從小區(qū)門口出發(fā)前500米完全重疊然后分道揚鑣一條去地鐵站一條去菜市場一條去公園。把所有點混在一起聚類前500米的點一定被劃成同一簇因為它們在空間上緊密抱團而后面的點按照目的地被切成兩三簇。最后你得到的結(jié)果名為“軌跡簇”實為“地理熱點”。如果你想分析的是“早上從小區(qū)出發(fā)的人分別去哪”這個結(jié)果當然可用但如果你想分析的是“三條完整軌跡的走法結(jié)構(gòu)”這種聚類方式從根上就錯了。1.3 采樣率不一致同名點對之間的“假距離”讓Kmeans無從下手就算你統(tǒng)一了軌跡長度還有一個隱蔽問題采樣率不一致。假設(shè)兩條軌跡都從A地到B地A用戶的車是10Hz高頻率上報B用戶的車是1Hz低頻上報。兩者重采樣到同樣50個點之后理論上看起來可比了但如果中間有轉(zhuǎn)彎高頻軌跡記錄下了細膩的弧線低頻軌跡則是一條折線對應點之間的空間距離會非常大。Kmeans算歐氏距離時這種誤差會被放大結(jié)果就是同一條路線被分成兩簇僅僅因為采樣率不同。我的判斷標準很簡單如果你的聚類目標是“整條軌跡的形態(tài)或行為”Kmeans就不能直接吃原始坐標點如果目標是“找出經(jīng)常被經(jīng)過的區(qū)域”那直接用坐標點聚類沒問題。想清楚這一點你才敢往下走。2. 軌跡聚類的三條技術(shù)路線選路線比寫代碼更重要前面說的都是“不能怎么做”接下來聊“該怎么做”。軌跡聚類沒有唯一標準答案我在實際項目里一般會根據(jù)數(shù)據(jù)條件和業(yè)務(wù)訴求在三條技術(shù)路線里選一條。這三條路線在工程里都非常常見選對了能省掉大量返工時間。2.1 路線A弧長重采樣加展平特征向量這是本文Matlab代碼采用的主路線也是我給大多數(shù)項目做“第一版結(jié)果”時的首選。核心思路分成三步對每條軌跡做弧長重采樣統(tǒng)一成N個坐標點把N個點的橫縱坐標拼接成一個1行2N列的向量作為這條軌跡的特征對這個特征矩陣跑標準Kmeans。這個方案最大的好處是簡單粗暴且可解釋每條軌跡的每個重采樣點都有固定含義特征向量里的第i個位置就是“整條軌跡在里程第i段處的坐標”。Kmeans聚類后簇心向量可以重新折疊成一條軌跡這就是這個簇的“平均軌跡”畫出來非常直觀。缺點是對軌跡的形變和錯位比較敏感兩條形態(tài)相似但起點終點偏移較大的軌跡可能因為末端坐標差過大被拆開。2.2 路線B距離矩陣加MDS或譜聚類如果你覺得軌跡形態(tài)差異很明顯起點終點都不固定用路線A效果不好那我建議用路線B。先算任意兩條軌跡之間的相似度距離得到一個距離矩陣D再通過多維縮放或者譜聚類把軌跡分組。這個路線不要求軌跡等長你可以用DTW距離、Hausdorff距離、最長公共子序列距離甚至自定義業(yè)務(wù)距離非常靈活。缺點是計算復雜度和存儲開銷都不小。60條軌跡的距離矩陣是60×60600條軌跡就是600×600如果每條軌跡的DTW計算都要幾十毫秒整體就是幾分鐘起步。所以這個路線更適合離線分析不適合在線實時聚類。2.3 路線C手工特征壓縮軌跡再跑Kmeans如果業(yè)務(wù)上面要求“每個簇必須講得出商業(yè)含義”那路線C往往是最明智的選擇。不要拿整條軌跡去聚類而是先從軌跡里提取出有業(yè)務(wù)含義的標量特征全程平均速度、最大瞬時速度、總里程、轉(zhuǎn)彎次數(shù)、平均航向角變化率、起終點直線距離與里程的比例等。把這些特征拼成一個向量再做Kmeans。這種做法的好處是聚類結(jié)果可以直接落到業(yè)務(wù)動作上。“簇1是速度高、轉(zhuǎn)彎少的長途直行軌跡簇2是速度低、轉(zhuǎn)彎多的短途巡游軌跡”這句話講給運營聽對方一下就懂了。缺點是特征怎么設(shè)計很考驗經(jīng)驗特征選得不好聚類結(jié)果等于隨機分組。2.4 三條路線怎么選一張表講清楚場景特征推薦路線理由軌跡長度差不遠采樣率經(jīng)過預處理路線A實現(xiàn)成本最低結(jié)果最穩(wěn)定軌跡形態(tài)差異大起點終點各不相同路線B距離度量更靈活能捕捉形變業(yè)務(wù)要求解釋性要直接給運營用路線C特征向量語義清晰可落行動數(shù)據(jù)量在幾百條以下離線分析路線A或B計算量可控Debug方便數(shù)據(jù)量上萬條在線聚類路線A特征固定后可直接用增量式Kmeans3. 預處理是命門弧長重采樣、數(shù)據(jù)清洗、坐標歸一化的細節(jié)很多人在Kmeans上報錯“矩陣維度不一致”就開始換框架、換語言其實問題出在前面處理原始軌跡時沒做好標準化的動作。預處理不是可有可無的步驟它對聚類結(jié)果的影響比Kmeans參數(shù)大一個數(shù)量級。3.1 弧長重采樣為什么必須按里程等間隔取點弧長重采樣的原理說白了就是按照軌跡的總長度把整條軌跡等分成若干段在等分點處取坐標。這樣做的好處是徹底消除采樣率和速度差異的影響。拿兩條從公司到家的軌跡舉例一條是早高峰堵車車速慢GPS點位密集另一條是半夜暢通車速快GPS點位稀疏。如果不重采樣點位密度直接反映車速Kmeans會誤把“堵車時段”當成一個聚類特征重采樣之后每條軌跡都用等量的點描述相同的地理距離比較的就只剩下“路線形態(tài)”了。Matlab里實現(xiàn)弧長重采樣很直接我自己常用的函數(shù)長這樣function trajN resampleTrajectory(traj, N) % traj: n x 2 的軌跡點列例如經(jīng)投影后的平面坐標 % N: 重采樣點數(shù)也就是統(tǒng)一后每條軌跡的點數(shù) % trajN: N x 2 的重采樣結(jié)果按全程弧長等間隔取點 if size(traj, 1) 2 error(軌跡點數(shù)太少無法重采樣); end segLen sqrt(sum(diff(traj, 1, 1).^2, 2)); cumLen [0; cumsum(segLen)]; totalLen cumLen(end); if totalLen 0 trajN repmat(traj(1,:), N, 1); return; end tq linspace(0, totalLen, N); % interp1 會對矩陣的每一列分別插值 trajN interp1(cumLen, traj, tq, linear); end關(guān)于重采樣點數(shù)N怎么選我的經(jīng)驗是看軌跡里最短的那條。最短軌跡如果有20個點那你重采樣到30~50個點是安全的信息不會丟太多如果最短軌跡只有5個點那這數(shù)據(jù)本身質(zhì)量就存疑硬重采樣到50只會把噪聲插值得漂漂亮亮沒意義。3.2 清洗跳變點和NaN先刪臟數(shù)據(jù)再談聚類GPS軌跡常見的臟數(shù)據(jù)有兩大類一類是NaN丟星另一類是跳變點。跳變點的典型表現(xiàn)是相鄰兩個采樣點之間距離超大比如1秒內(nèi)“瞬移”了500米這通常是GPS漂移造成的。清洗跳變點的方法也比較樸素先計算每相鄰兩點的速度設(shè)定一個閾值比如超過整條軌跡速度中位數(shù)的5倍就判定為跳變把跳變點刪除再用前后有效點線性插值補上。Matlab里判斷跳變的核心代碼段我貼出來dt diff(t); % t 是時間序列單位秒 segSpeed sqrt(sum(diff(xy, 1, 1).^2, 2)) ./ dt; medSpeed median(segSpeed); outlierIdx find(segSpeed 5 * medSpeed); % 設(shè)定5倍中位數(shù)為跳變閾值注意不要直接用絕對速度閾值因為不同場景的速度分布差太遠了。步行軌跡的中位數(shù)速度只有1.5m/s左右機動車軌跡的中位數(shù)速度輕松上10m/s用絕對閾值會誤殺。3.3 坐標歸一化比你想的更微妙如果你已經(jīng)在軌跡數(shù)據(jù)上做了重采樣那么最后的歸一化步驟直接影響Kmeans的聚類形狀。這里有個細節(jié)經(jīng)常被忽略到底是對每一列單獨zscore還是對全部坐標做統(tǒng)一縮放我的建議是如果你要保留軌跡的形狀比例就做統(tǒng)一縮放不要對x坐標和y坐標分別標準化。分別標準化等于把橫軸和縱軸的尺度強行拉成一樣圓形軌跡會被壓成橢圓直角三角形會被壓成斜邊變短的畸變?nèi)切巍Ｕ_的做法是先求出所有軌跡在x和y方向上的整體標準差然后用同一個尺度因子縮放。實際操作里我更多是直接對展平后的特征矩陣做一次統(tǒng)一標準化X reshape(trajRes, numTraj, N * 2); scaleFactor std(X(:)); X X ./ scaleFactor;這樣每個特征維度的相對比值還保留著原始的物理意義同時又讓量級進入了Kmeans友好的范圍。4. 相似度度量RMSD、DTW與Hausdorff的取舍邏輯Kmeans的實質(zhì)是基于距離的迭代優(yōu)化所以距離度量才是那個真正決定聚類結(jié)果的東西。很多人的代碼跑得很熟練但從不回頭看一眼自己用的是哪種距離出了奇怪結(jié)果也不知道從哪查起。4.1 RMSD路線A背后的默認度量在重采樣之后的軌跡空間里最自然的距離是逐點距離的平均值也就是RMSDRoot Mean Square Deviation。兩條重采樣軌跡都是N×2的矩陣相減之后求每個點的歐氏距離然后取平均。d sqrt(sum((trajA - trajB).^2, 2)); dRMS mean(d);RMSD的好處是計算快而且對小幅噪聲有天然平滑效果。但它對時間軸錯位非常敏感如果一條軌跡在某個彎道比另一條早了10米開始轉(zhuǎn)彎逐點對應時轉(zhuǎn)彎前的直線段會被強行錯位比較距離值虛高。這也是為什么路線A對“形變”不友好。4.2 DTW路線B最常用的距離代碼不長但邏輯要懂DTW動態(tài)時間規(guī)整解決的就是RMSD對時間軸錯位的毛病。它的思想是允許軌跡點按順序“錯位對齊”代價最小的對齊方式就是最優(yōu)匹配。比如兩條軌跡同樣是繞個S彎一條彎得早、一條彎得晚DTW能夠正確地把兩邊的彎心匹配上距離自然就小了。Matlab實現(xiàn)版本很多我自己用的簡潔版本如下function d dtwDist(tr1, tr2) % tr1, tr2: n1 x 2, n2 x 2 的兩條軌跡點列 n1 size(tr1, 1); n2 size(tr2, 1); D pdist2(tr1, tr2); % 點對點距離矩陣 C inf(n1 1, n2 1); C(1,1) 0; for i 1:n1 for j 1:n2 % 只能向前對齊不能回頭 C(i1,j1) D(i,j) min([C(i,j1), C(i1,j), C(i,j)]); end end d C(n11, n21); end這個雙重循環(huán)實現(xiàn)直觀但壞消息是復雜度是O(n1×n2)軌跡一長、樣本一多就跑不動。實際項目里我一般只在軌跡數(shù)量不超過200條時用DTW全量距離矩陣。如果想提速可以限制扭曲窗口寬度比如只允許點對在前后20%長度范圍內(nèi)對齊效果幾乎不變速度能快好幾倍。4.3 Hausdorff距離和Frechet距離一句話講清楚差異在路線B里還可以用Hausdorff距離或Frechet距離。Hausdorff距離的定義是“軌跡A上每個點到軌跡B的最近距離的最大值”它的特點是極其關(guān)注兩條軌跡的最大偏離程度。這個性質(zhì)很雙刃它擅長找出“某一段明顯背離”的軌跡但也很容易被單個漂移點帶偏。Frechet距離則更接近“人在遛狗狗繩拉直時最短能多短”的模型它同時約束了點位的順序性和連續(xù)性比Hausdorff更符合軌跡形態(tài)比較的直覺。但Frechet的計算比DTW更復雜Matlab沒有內(nèi)置函數(shù)需要自己實現(xiàn)離散版本工程成本高一般我不推薦第一個版本就上它。4.4 實際該用哪個結(jié)合你的數(shù)據(jù)特性拍板如果你的軌跡已經(jīng)做了弧長重采樣并且采樣率差異不大直接用RMSD就行Kmeans的收斂速度快結(jié)果也好解釋。如果軌跡在空間上有明顯的“先經(jīng)過A再經(jīng)過B”的順序結(jié)構(gòu)但時間節(jié)奏不同用DTW。如果只想粗篩“有沒有某一段軌跡嚴重偏離”用Hausdorff。如果業(yè)務(wù)上把軌跡當成“一條繩”要求整條路徑的形狀逼真度再考慮Frechet。5. 完整Matlab實現(xiàn)從模擬軌跡到Kmeans聚類的一站式腳本接下來是這篇文章的重頭戲一份可以直接運行的Matlab腳本。為了讓你不依賴外部數(shù)據(jù)就能看到完整效果我先生成三組形態(tài)不同的模擬軌跡然后做弧長重采樣、特征展平、Kmeans聚類和可視化。你把自己真實的軌跡數(shù)據(jù)套進對應的步驟就行。5.1 主腳本全流程一口氣跑通%% 基于Kmeans的軌跡聚類完整示例 clear; clc; close all; rng(42); %% Step 1: 模擬生成三種形態(tài)的軌跡 numTraj 90; % 總軌跡條數(shù) trajCell cell(numTraj, 1); for i 1:numTraj if i 30 % 第一類平緩直線型 nPts 18 randi(10); t linspace(0, 1, nPts); x 0.8*t 0.04*randn(nPts,1) 0.1*randn; y -0.3*t 0.04*randn(nPts,1) 0.1*randn; elseif i 60 % 第二類上凸曲線型 nPts 16 randi(12); t linspace(0, 1, nPts); x t.^2 0.05*randn(nPts,1); y sin(1.5*pi*t) 0.06*randn(nPts,1) 0.1*randn; else % 第三類先直后折的折線型 nPts 24; t linspace(0, 1, nPts); x min(0.6*t, 0.6) 0.03*randn(nPts,1); y 0.8*t 0.03*randn(nPts,1); y(t 0.6) 0.48 1.2*(t(t0.6)-0.6) 0.03*randn(sum(t0.6),1); end trajCell{i} [x, y]; end %% Step 2: 弧長重采樣到統(tǒng)一點數(shù) N 50; % 統(tǒng)一后的軌跡點數(shù) numTraj length(trajCell); trajRes zeros(numTraj, N, 2); for i 1:numTraj trajRes(i,:,:) resampleTrajectory(trajCell{i}, N); end %% Step 3: 展平成特征向量并做統(tǒng)一縮放 X reshape(trajRes, numTraj, N*2); scaleFactor std(X(:)); X X ./ scaleFactor; %% Step 4: 用輪廓系數(shù)掃描K值K從2到6 Krange 2:6; silScores zeros(1, length(Krange)); for k 1:length(Krange) idxTmp kmeans(X, Krange(k), Replicates, 10); silTmp silhouette(X, idxTmp); silScores(k) mean(silTmp); end [~, bestKPos] max(silScores); K Krange(bestKPos); %% Step 5: 用最優(yōu)K值跑最終Kmeans [idx, C] kmeans(X, K, Replicates, 20); %% Step 6: 分別計算每個簇的平均軌跡用于可視化簇中心 avgTrajCell cell(K, 1); for k 1:K members find(idx k); avgTrajCell{k} squeeze(mean(trajRes(members, :, :), 1)); end %% Step 7: 可視化 figure(Position, [100 100 1200 420]); colors lines(K); subplot(1,3,1); hold on; grid on; for i 1:numTraj traj_i squeeze(trajRes(i,:,:)); plot(traj_i(:,1), traj_i(:,2), Color, [0.65 0.65 0.65]); end title(重采樣后的所有軌跡); axis equal; subplot(1,3,2); hold on; grid on; for i 1:numTraj traj_i squeeze(trajRes(i,:,:)); plot(traj_i(:,1), traj_i(:,2), Color, colors(idx(i), :), LineWidth, 1.0); end title([Kmeans聚類結(jié)果K, num2str(K)]); axis equal; subplot(1,3,3); hold on; grid on; for k 1:K plot(avgTrajCell{k}(:,1), avgTrajCell{k}(:,2), ... Color, colors(k,:), LineWidth, 2.5); end title(每個簇的平均軌跡); axis equal;5.2 一個容易被可視化帶歪的細節(jié)簇中心來自標準化空間第一次用Kmeans做軌跡聚類的人很容易直接把kmeans返回的簇心C拿來畫“中心軌跡”畫出來的圖往往亂七八糟。原因在于我們喂給kmeans的特征矩陣X是經(jīng)過統(tǒng)一縮放的簇心C也是縮放空間里的坐標直接reshape回軌跡形狀時橫縱坐標已經(jīng)脫離了原始數(shù)據(jù)的物理尺度。所以在Step 6里我沒有用C去畫圖而是重新取出每個簇內(nèi)的原始重采樣軌跡按簇內(nèi)成員做平均。這樣得到的平均軌跡才是“這個簇的代表性路線”畫出來和原始坐標對得上業(yè)務(wù)人員看了也不會懵。這個細節(jié)我強烈建議你保留因為它直接影響你給同事或客戶匯報時的可信度。5.3 拿到聚類結(jié)果之后至少看一眼簇內(nèi)的平均軌跡代碼跑通只是開始。我的習慣是每次聚類結(jié)束第一件事不是看輪廓系數(shù)而是把每個簇的平均軌跡和兩三條典型單條軌跡疊在一張圖上快速確認“這個簇的代表性軌跡是否符合直覺”。如果平均軌跡雜亂無章說明K值選大了或者預處理出了問題如果平均軌跡平滑清晰才值得繼續(xù)往下做業(yè)務(wù)分析。6. K值選擇與聚類質(zhì)量評估輪廓系數(shù)、穩(wěn)定性與業(yè)務(wù)校驗Kmeans有個繞不開的宿命K要你自己給。很多教程講到這里就扔給你一個輪廓系數(shù)讓人背公式但我更愿意把K值選擇理解成“聚類結(jié)果的復現(xiàn)性校驗”。因為軌跡數(shù)據(jù)往往沒有標準答案同一個K在不同初始點下跑出來的穩(wěn)定程度才是更靠譜的評估指標。6.1 輪廓系數(shù)Matlab一行搞定但別只取一個值輪廓系數(shù)的計算方法你自己寫也不難對每個樣本算它到同簇其他樣本的平均距離a再算它到最近其他簇所有樣本的平均距離b輪廓系數(shù)就是(b-a)/max(a,b)。取值范圍從-1到1越接近1說明樣本離自己簇越近、離其他簇越遠。Matlab里直接調(diào)用內(nèi)置函數(shù)sil silhouette(X, idx); meanSil mean(sil);但我不建議只掃描一次就拍板K。腳本里我循環(huán)了K2到6然后取了平均輪廓系數(shù)最大的K。這里有個隱形問題輪廓系數(shù)容易偏袒K小的情況因為簇少的時候簇間邊界天然更清晰。所以我通常會把輪廓系數(shù)排名前兩三名都跑一遍再對比可視化結(jié)果選那個業(yè)務(wù)上最有意義同時輪廓系數(shù)也不差的K。6.2 肘部法則看簇內(nèi)距離平方和隨K的變化肘部法則是另一種常見的K值選擇方法看的是簇內(nèi)距離平方和WSS隨K增大的下降趨勢。下降速度驟減的那個拐點就是“肘部”。Matlab計算WSS可以手動實現(xiàn)for k Krange idxTmp kmeans(X, k, Replicates, 10); wss(k-1) 0; for kk 1:k members X(idxTmp kk, :); center mean(members, 1); wss(k-1) wss(k-1) sum(sum((members - center).^2)); end end plot(Krange, wss, o-);不過說實話真實軌跡數(shù)據(jù)里WSS肘部經(jīng)常不是那么明顯曲線平滑下降很難說哪個點才是肘。所以我的策略是肘部法用來圈定候選范圍輪廓系數(shù)用來進一步篩選最后用穩(wěn)定性拍板。6.3 穩(wěn)定性校驗同一K在不同初始化下是否總是給出相似分組這個指標很多教程不提但在業(yè)務(wù)項目里非常實用對同一個K用不同的隨機種子跑10次Kmeans然后比較10次結(jié)果的相似度。相似度高的說明這個K和這組數(shù)據(jù)是“匹配”的相似度低說明數(shù)據(jù)在這個K下本來就分不清任何結(jié)論都不可靠。Matlab里實現(xiàn)也不復雜用RandStream控制隨機種子跑10次聚類之后逐樣本比較標簽的一致性。如果標簽之間的平均互信息AMI低于某個閾值我就直接建議用戶換距離度量或者換路線而不是繼續(xù)調(diào)K。6.4 別讓數(shù)學指標凌駕于業(yè)務(wù)之上做軌跡聚類到最后我極少遇到“數(shù)學上的最優(yōu)K”和“業(yè)務(wù)上的最優(yōu)K”完全一致的情況。有一次共享單車項目的數(shù)據(jù)輪廓系數(shù)最高指向K6但運營那邊實際只有三種可執(zhí)行的調(diào)度策略K6意味著每種策略要被拆成兩簇操作上根本無法落地。后來我折中選了K4在多簇之間加了細分標簽輪廓系數(shù)從0.55掉到0.47但每個簇都能對應到明確的運營動作。這件事給我的教訓是K值選擇歸根結(jié)底是服務(wù)業(yè)務(wù)目標的決策數(shù)學指標是工具不是判決書。7. 真實軌跡數(shù)據(jù)上的踩坑記錄采樣率、停留點與離群點的處理最后這部分是這些年在真實項目里積累下來的一些經(jīng)驗。模擬數(shù)據(jù)再怎么完美到了真實軌跡場景該遇到的坑一個都少不了。我把最典型的幾個記錄下來希望你不用重新踩一遍。7.1 采樣率差異的威力以出租車GPS數(shù)據(jù)為例出租車GPS數(shù)據(jù)的采樣間隔并不是恒定的空載省電模式下可能是30秒一條載客接單后變成5秒一條。兩條從機場到市中心的軌跡如果一條是空載狀態(tài)錄的一條是載客狀態(tài)錄的重采樣之前算距離數(shù)值可以差到幾公里?；￠L重采樣能解決點位數(shù)量不均的問題但沒法徹底解決“重要拐彎處點數(shù)少”的問題。我通常會額外做一步先按時間插值到1秒間隔再做弧長重采樣。這樣等于先補密再統(tǒng)一拐彎處的形態(tài)特征保留得更好。7.2 停留點會扭曲軌跡形狀必須識別并剁掉軌跡里最坑人的是停留行為。用戶在某地停了20分鐘GPS以1Hz頻率一直在上報這段軌跡在空間上表現(xiàn)為一個點團。如果你把整條軌跡拿去重采樣這段點團會占掉重采樣點里很大比例于是這條軌跡的“形狀”被這個停留點完全主導。聚類時不管它真正的行駛路線是什么樣子都會被分到“停留時間長的軌跡”那一簇去。我的處理方式是先算每個相鄰點對的速度把速度低于0.5m/s的連續(xù)片段標記為停留段然后將停留段內(nèi)的點做降采樣比如每20個點保留1個或者干脆把停留段的中心點作為單點保留。具體保留策略要看業(yè)務(wù)如果分析的是路徑規(guī)劃停留點直接刪掉如果分析的是出行行為停留時長本身反而是一個特征那就把它提取出來不要混在坐標里。7.3 離群點Kmeans最怕那種“又長又怪”的軌跡離群軌跡對Kmeans的影響比離群點對大得多。比如100條正常軌跡里混進來一條把整個城市的對角線都跑了一遍的軌跡它的長度和形狀都極端重采樣之后這個樣本在特征空間里距離其他樣本極遠。這會導致兩個后果要么它單獨成一簇K被它浪費掉一簇要么它強行拽動某個簇心把正常軌跡也帶偏。我處理離群軌跡的順序是先用簡單的長度和里程閾值過濾掉明顯異常的軌跡再用DBSCAN在展平特征空間上做一次粗聚類把落在任何簇外或者簇很稀疏的樣本標記為離群。最后才把這些干凈樣本送入Kmeans。這樣不僅結(jié)果穩(wěn)還能在報告里多出一個“離群軌跡”類別很多情況下這個類別反而能發(fā)現(xiàn)異常駕駛或設(shè)備故障比正常聚類結(jié)果更有業(yè)務(wù)價值。7.4 在線跑Kmeans時的坑中心漂移與增量更新最后提醒一下想做在線軌跡聚類的朋友。Kmeans的原始版本是批處理的每來一批新軌跡就要重新跑一遍全量數(shù)據(jù)隨著數(shù)據(jù)量增長越來越慢。實際工程中我多數(shù)時候改用MiniBatchKMeans或者在線式增量更新新樣本進來后先算它到現(xiàn)有簇心的距離歸入最近簇同時按一定學習率更新簇心。Matlab里沒有內(nèi)置MiniBatchKMeans但你可以自己寫一個幾十行的更新循環(huán)核心就是那個簇心更新公式。如果你的軌跡流是實時上報的這個方向會比反復調(diào)用全量kmeans靠譜得多。……關(guān)于代碼最后再補充一點如果你自己只有經(jīng)緯度坐標建議先把經(jīng)緯度轉(zhuǎn)換成平面坐標再跑聚類。Matlab的Mapping Toolbox有deg2km或geodetic2enu這類函數(shù)沒有工具箱的話用等距投影的近似公式也夠用。這個步驟不做你在真實地理數(shù)據(jù)上計算歐氏距離誤差在低緯度地區(qū)還勉強能接受到了高緯度地區(qū)單位經(jīng)度對應的地面距離變化非常大聚類結(jié)果會失真到?jīng)]法看。我個人的習慣是拿到軌跡后的第一件事永遠是畫圖先把所有軌跡疊一張圖看一遍再決定用什么距離、什么路線、什么K。這個習慣幫我避開了至少一半的無效調(diào)試希望你也能用上。