手機信令的時空密度軌跡點識別算法

陳 略，熊 宸，蔡 銘

（中山大學智能工程學院廣東省智能交通系統重點實驗室，廣州 510006）

0 概述

近年來，手機信令作為一種持續時間長、覆蓋范圍廣的新型大數據被廣泛用于交通出行行為和規律的研究。手機信令數據量大、采樣頻率不均、定位精度低以及基站振蕩等問題，給其研究帶來諸多難題。要實現手機信令在現實中的廣泛應用，針對手機信令數據的基礎性研究尤為關鍵。手機信令預處理與軌跡點分析作為手機信令研究的第一步，其分析結果是研究用戶出行行為的基礎，對研究用戶出行行為與規律具有重要意義。

本文通過闡述手機信令特性，總結現有對數據預處理與停留點識別的研究成果，提出一種手機信令時空密度軌跡點識別算法。采用網絡軌跡時空聯結消除手機信令的空間不確定性，融合軌跡的曲折性、移動時間和停留時間重新定義網格簇內軌跡點時空移動能力，以網格簇時空密度為特征識別停留簇，同時使用基于簇的時空關系漸進聚類算法合并時空鄰近的停留簇，并將其識別為停留區域。

1 相關工作

1.1 手機信令的特性

手機信令是一種時間序列數據，由基站的經緯度和時間戳構成。手機信令的產生有主動產生和被動產生兩種方式。主動產生方式主要包括用數據流量上網、發短信、打電話等主動行為觸發的基站響應；被動產生方式主要包括收短信、接電話等被動行為觸發的基站響應。由于手機數據的產生無規律性，因此手機信令存在采樣頻率不均［1］、軌跡點時間間隔不規則以及軌跡點密度分布不均等問題。如果以0.5 h 為單位，將一天24 h 劃分為48 個區間［2］來統計手機信令軌跡數據的區間數，則一天中未達到16個區間分布的軌跡數據視為稀疏數據。由于手機信令采樣頻率不均，在手機主動上網或移動的過程中，手機信令的采樣間隔最小為1 s，并會連續觸發同一個位置的基站，造成未經預處理的手機信令數據量大且信息冗余，不適合直接進行軌跡點識別。

手機信令相較GPS、卡口數據具有更大的空間不確定性，該特性由手機信令的定位精度和工作機制所決定。用戶位置一般用其所在基站的位置表示，由于基站覆蓋范圍從幾百米到幾千米不等［3］，因此手機信令精度較差。同時，受手機信令工作機制的影響，用戶所在位置實際連接的基站可能并不是與其最接近的基站，由于某些區域的基站密度較大，某個地理位置可能同時被幾個基站所覆蓋，因此在該區域手機信號常隨基站信號強度變化而不斷切換基站，出現乒乓切換現象。一般而言，基站有負荷優化調節機制，當鄰近基站用戶負荷過大時，會自動切換到更遠但用戶負荷較少的基站，從而出現信號漂移［4］現象。此外，地形起伏、高樓遮擋等因素還會造成異常切換，從而引起手機在兩個或多個基站之間相互切換的乒乓效應或信號漂移現象。據文獻［5］統計，出現乒乓切換和信號漂移的數據量約占總數據量的30%。

乒乓切換和信號漂移并不代表用戶存在真實的移動，如果不及時對其進行處理，會因高速度和遠距離的切換或漂移影響停留點聚類，導致產生許多虛假出行數據。因此，針對乒乓切換和信息漂移造成的空間不確定性，研究人員嘗試用速度或距離清洗規則清除含有該特征的數據，但由于乒乓切換與信息漂移規律尚不明確，清洗規則依據不足，且清除效果未知，因此通過清洗來消除空間不確定性還有待研究。本文采用不同于上述方法的思路，不以清除所有乒乓切換和漂移數據為目的，而將軌跡點震蕩特征作為時空鄰近一簇軌跡點徘徊與潛在停留的標志，并將具有這種特征的軌跡點通過時空關系聯結為軌跡簇。通過將震蕩的軌跡點同化為一個整體，忽略簇中軌跡點的經緯度和時間戳，并在算法中以簇為聚類單元，消除乒乓切換和漂移數據產生的空間不確定性。

手機信令具有空間不確定性，其除了是軌跡數據之外，還是一種典型的時空數據［6-8］。這類軌跡在聚類時需要同時考慮時間和空間兩個維度，而目前缺乏以時間和空間兩個維度度量軌跡點停留屬性強度的指標，大部分算法仍采用時空分離的聚類算法，按照先時間后空間或者先空間后時間的順序進行聚類。本文結合軌跡點簇內部以及前后軌跡簇間的時間、距離和曲折性定義軌跡簇的移動能力，并根據軌跡簇移動能力及其與前后軌跡簇的時空距離提出軌跡簇時空密度函數MAST 來度量軌跡簇停留屬性的強度，并為停留點識別提供具體的參照指標。

1.2 手機信令預處理

手機信令預處理包括對時間連續的同位置點進行合并以減少軌跡量，以及對原始數據中大量乒乓切換數據和漂移數據進行處理，以消除空間不確定性。目前對乒乓切換數據和漂移數據的處理［4，9］有分層清洗［10-11］和多階段去噪［1］兩種方式。

分層清洗是先將噪聲分為乒乓切換數據、漂移數據等，再針對各類噪聲的特征制定不同規則進行清洗。乒乓切換數據主要有速度規則［12］和模式規則［5，13］兩種清洗規則。速度規則為：若基站之間存在ABA 切換形式，且基站間速度超過一定閾值，則認為是乒乓切換，可按照一定規則清除。在模式規則方面，文獻［11］通過定義乒乓切換形式，設定在基站間切換的模式，若檢測到該模式則進行清除。文獻［14］中結合速度規則和模式規則的乒乓切換清洗也較常見。對于漂移數據，通常認為其在短時間內與前后兩個軌跡點相距較遠，一般采用設定中間軌跡點與前后軌跡點間的速度來檢測漂移數據［10］，如果速度超過閾值，則認為是漂移數據；或者檢測漂移點與前后軌跡點形成的角度，通過設置角度閾值來清洗漂移數據［11］。

多階段去噪是將乒乓切換數據和漂移數據導致的偏差視為由移動通信機制和采集系統故障引起的系統誤差，這種誤差難以用數學工具來描述，即無法用任何隨機分布形式來模擬和表達，但可將其在二維上表現為離群點，采用離群點檢測算法進行處理［1］。與制定多種檢測規則來清洗數據的分層清洗不同，多階段去噪是將乒乓切換數據與漂移數據統一視為離群點進行處理。實際上，由于乒乓切換和漂移效應的形式具有多樣性，且目前乒乓切換數據與漂移數據的內在規律尚不明確，因此無法評價清洗的有效性。文獻［15］借鑒貪婪算法思想根據位置出現頻率提出時空貪婪同化算法，尋找短時間內該位置點之間夾雜的其他位置點并形成集合，以實現位置點同化，從而消除空間不確定性，錨固長時間停留點的位置。但是僅按照頻率自上而下合并，缺乏對于不同時間段的同一位置點（手機信令時間序列性）的考慮，容易將某些經常移動的點識別為停留點從而打亂軌跡序列。

消除空間不確定性在聚類中表現為受到噪聲點的影響較小。由各種消除空間不確定性的算法可知：清洗思想以軌跡點為研究單元，通過劃分時間閾值和距離閾值進行聚類，易被乒乓切換噪聲點及漂移噪聲點割裂為幾個小簇；同化思想是將有震蕩和徘徊等潛在停留特征的軌跡點同化為簇，并以簇為研究單元在軌跡簇之間進行聚類，受噪聲影響較少。

由上述可知，同化算法更適合應用于噪聲多、精度低以及復雜度大的手機信令數據。針對現有同化算法未考慮手機信令時間序列性，將高頻率位置點間夾雜的低頻率位置點歸為一簇，并將所有時刻的該位置點歸于其軌跡簇而忽視軌跡復雜性等問題，本文根據乒乓切換和數據漂移共有的震蕩特征及交通出行時間的定義，設計在一定時間內空間鄰近位置點及其夾雜的位置點互相環扣進行交織聯結的時空聯結流程，同時兼顧軌跡的時序性與時空緊密性。

1.3 停留點識別

停留點是一次出行的起訖點，其通常被定義為有目的的活動地點，且該地點的停留時間超過5 min。由于停留時間需要連續的空間累積，因此停留點的識別算法包括基于距離、時間、形狀和密度的聚類算法。

基于距離的聚類算法［16-17］基于手機信令采樣頻率不固定的特征，通過相鄰軌跡點間的距離來反映用戶出行情況，以代替用相鄰軌跡點間的時間差作為停留時間來判斷用戶是否停留。因此，可根據前一個軌跡點的移動或停留狀態，以及前一個軌跡點與當前軌跡點的距離與時間差判斷當前軌跡點的狀態。但這種基于距離的算法易受到漂移點和前一個軌跡點狀態判斷錯誤的影響。此外，距離閾值和時間閾值的設定對聚類結果也有較大影響。

基于時間的聚類算法［18］沿時間軸將距離小于閾值的軌跡點聚類成簇，并通過判斷簇的停留時間來識別是否為停留區域。基于時間的聚類算法易受漂移點和震蕩點影響，從而將停留區域切分為多個小簇，且由于沿時間軸進行聚類，會將真正引起簇內距離過大的時間軸靠前的軌跡點納入簇中，造成簇的劃分錯誤。因此，時間閾值和距離閾值的選取在很大程度上會影響聚類劃分結果。

基于軌跡形狀［11］用定義的最小覆蓋圓沿時間軸將軌跡點劃分到簇中，當軌跡點超過最小覆蓋圓時判斷簇的停留時間，若大于停留時間閾值則判斷為停留簇，若小于停留時間閾值則排除簇中時間軸最靠前的軌跡點，并從簇中第二個軌跡點開始循環。基于形狀的聚類算法不受聚類順序的影響，較基于時間聚類的算法能更有效地劃分距離近的簇，但容易受到噪聲影響而分為多個小簇，且該算法受停留時間閾值的影響較大。

DBSCAN 等基于密度的聚類算法根據手機信令的時間序列性，將Eps 空間鄰域拓展到時間維度［10，19-20］。時空密度被定義為時空閾值范圍內的軌跡點個數，對于采樣頻率不均的手機信令數據而言，不能將軌跡個數作為時空密度的度量標準，加上密度聚類參數選取困難，當采樣頻率不均時，聚類質量較差。

目前軌跡聚類算法存在以下問題：

1）對停留點構建指標的判斷未結合時空兩方面屬性，只能定性判斷軌跡點停留與否，無法定量反映軌跡點停留屬性的大小。

2）僅將軌跡作為線性曲線處理，未考慮軌跡形狀、曲折性等其他復雜的特性。

3）對于時間、距離、個數等聚類算法的閾值采取統一方式，未兼顧不同用戶的出行特征。

針對上述問題，本文提出解決方案如下：

1）定義了融合時空特性的時空密度函數MAST指標，其用來度量軌跡簇中軌跡點的時空密度，并定量反映軌跡點停留屬性的強度。

2）根據軌跡簇的曲折性以及移動時間和停留時間比例定義簇的移動能力，并作為軌跡特性融合到時空密度函數計算中。

3）不同于聚類算法中的閾值設置方法，本文根據時空密度函數MAST 的概率密度分布特征，自動劃分符合個人移動特征的停留點密度閾值，避免了手動設置閾值不適用于全局數據的問題。

2 本文算法

本文提出基于時空聯結移動能力的軌跡點識別算法，其流程如圖1 所示。

圖1 本文算法流程Fig.1 Procedure of the proposed algorithm

2.1 軌跡點網格化

為便于統一算法評估尺度與簡化計算量，將研究區域網格化，并將基站的經緯度歸入網格中，以網格序號代替。例如，將第i條記錄軌跡點基站的經緯度(i，lngi，lati，ti)替換為網格坐標(i，xi，yi，ti)。

將研究區域網格化后，對連續的同網格軌跡點記錄進行合并，只保留第一條和最后一條。例如，序列{(i，x，y，ti)，(i+1，x，y，ti+1)，…，(j，x，y，tj)}簡化為{(i，x，y，ti)，(j，x，y，tj)}。連續的網格軌跡點記錄合并后成為前后兩條相同網格的時間戳記錄，由此可以明確在此兩個時間戳及其之間時段內該網格的停留時間，并將該網格視為停留網格。由于部分網格無連續軌跡點記錄，只有單個網格的時間戳，因此將此類網格視為單網格，無法知道該網格的停留時間。

2.2 網格軌跡時空聯結

本文時空聯結的思想是基于手機信令的振蕩噪聲及其停留與徘徊的特點：從起點切換到某地再回到起點附近的過程如果發生在短時間內，那么這段時間用戶實際并未出行，僅在起點和切換點之間的某處靜止或小范圍地運動。由各城市居民出行的調查結果可知，一次出行的耗時通常大于5 min，其中在出行目的地的耗時超過5 min。由于從某地點出發再返回該地點附近的過程涉及一個活動和兩次出行［15］，因此該過程的最小時間間隔為15 min。

算法1網格軌跡時空聯結算法

上述算法具體步驟如下：

1）假設網格化后軌跡數據集為D，cluster 為標識的網格記錄集合，del_cluster 為待循環刪除的網格記錄集合，Cm 為所有網格簇集合，當前軌跡點記錄(i，xi，yi，ti)∈cluster。

2）按時間序列輸入當前軌跡點記錄(i，xi，yi，ti)，若下一條記錄(i+1，xi+1，yi+1，ti+1)為相鄰網格，則xi+1∈[xi-1，xi+1]且yi+1∈[yi-1，yi+1]。若(i+1，xi+1，yi+1，ti+1)不在cluster 中，則將下一條記錄(i+1，xi+1，yi+1，ti+1)加入簇和del_cluster中，當前記錄變為下一條記錄(i+1，xi+1，yi+1，ti+1)，即i=i+1，然后再從步驟2 開始；否則轉到步驟4。若下一條記錄所在網格與當前網格不相鄰，則轉到步驟3。

3）檢測ti＜t＜ti+15 min 時間范圍內是否有與當前軌跡點網格相鄰的記錄，即xi+1∈[xi-1，xi+1]且yi+1∈[yi-1，yi+1]，將與當前軌跡點網格相鄰的記錄集合記為cf。如果該時間段內存在相鄰網格記錄，即length（cf）>0，則檢測相鄰網格記錄是否存在不屬于cluster 的軌跡記錄，并將其記為ss，如果存在該軌跡記錄，即length（ss）>0，則將這部分不屬于cluster的相鄰網格軌跡記錄及其時間范圍內所有不屬于cluster 的軌跡記錄全部加入cluster 和del_cluster，當前記錄轉為cluster 中最后一條記錄，即i=cluster［-1］.order number。若不存在不屬于cluster 的軌跡記錄，則轉到步驟4，執行算法2；若ti＜t＜ti+15 min 時段內沒有與當前軌跡點網格相鄰的記錄，即length（cf）=0，則轉到步驟4，執行算法2。

4）刪除del_cluster 中所有與當前軌跡網格序號相同（x=xi且y=yi）的記錄，若del_cluster 仍不為空，則按時間排序，當前軌跡記錄跳轉到del_cluster 中時間排序最后的軌跡記錄，即i=del_cluster［-1］.order number（見算法2 中步驟1～步驟3）。若刪除當前軌跡記錄所在網格后del_cluster 為空，則轉向步驟5。

5）將cluster 中網格軌跡記錄標記為潛在的停留簇，并將cluster 加入Cm，當前記錄跳轉到cluster 中最后一條記錄的下一條，即i=cluster［-1］.order number+1（見算法2 中步驟5～步驟6）。

算法2刪除與循環聯結算法

算法2 是網格軌跡時空聯結算法的子算法，其執行的是算法1 中在當前軌跡點不能再向后聯結軌跡點的情況下，在簇軌跡點中從后向前尋找能繼續沿時間序列向后聯結的軌跡點作為當前軌跡點的過程。若簇內沒有可以繼續向后聯結的軌跡點，則結束該簇聯結流程，并將簇加入到Cm 中。算法2 的具體步驟為：在當前軌跡點與下一個軌跡點的網格位置不相鄰，且當前軌跡點時間戳后的15 min 內無相鄰網格可以聯結情況下，刪除del_cluster 簇中全部當前網格序號的軌跡點，并在簇del_cluster 剩余的網格軌跡點中選擇時間戳最后的軌跡網格作為當前網格，繼續對該簇循環執行算法1 的聯結流程。如果刪除后del_cluster 為空，無法再向后聯結，則結束該簇的聯結流程，將簇cluster 添加到Cm 集合中，向后聯結下一簇。

在時空聯結的簇中，停留網格記錄需成對出現。通過時空聯結的循環流程，可按照時間序列將一定時段和地理范圍內停留或徘徊交織的軌跡點簇聯結為網格簇。時空聯結流程偽代碼見算法1 和算法2。經過時空聯結，網格主要有單網格、停留雙網格以及網格簇3 種類型。

2.3 網格時空密度

本文參考文獻［21-22］中密度的定義提出一種新的時空密度函數MAST，根據軌跡段的移動曲線特征并結合距離和時間分布情況，定義網格簇的移動能力和網格簇的時空密度，以體現出信令的時空特性與時間序列性。

圖2為兩種軌跡的移動情形。可以看出在同一個起始點和目的地的軌跡段中，移動的軌跡偏向進行直線運動（見軌跡1），而停留或徘徊的軌跡更偏向進行曲線運動（見軌跡2）。文獻［21］根據圖2 提出移動能力的概念，基于此，下文對直線距離和曲線距離進行定義。

圖2 兩種軌跡移動情形Fig.2 Two cases of track movement

定義1（直線距離）存在n-m+1 個點的軌跡段Traji={pm，pm+1，…，pn}，Traji的直線距離為軌跡段起點pm和終點pn之間的直線距離，其表達式如下：

定義2（曲線距離）存在n-m+1 個點的軌跡段Traji={pm，pm+1，…，pn}，Traji的曲線距離為各子軌跡段的距離之和，其表達式如下：

移動能力被定義為直線距離和曲線距離的比值［21］，其表達式如下：

文獻［21］對MA 的定義僅考慮整個軌跡段的曲折性，而長時間在職住地或者基站覆蓋范圍較大的區域停留時，與其他基站的切換較少但停留時間較長，如果只考慮曲折性，那么在停留區域位置不變的點的移動能力為1，其移動能力反而會異常大。因此，本文綜合考慮軌跡的曲折性、移動時間和停留時間后，對網格時空移動能力進行重新定義。文獻［21］針對核心點通過該點所在軌跡段的MA 及其周圍網絡接入點（Nap）對其產生的影響來定義，但其利用周圍軌跡點個數劃定移動能力的計算范圍不適用于采樣頻率不均的手機信令數據。文獻［22］基于軌跡點鄰域半徑R內包含所有軌跡段的最小值對軌跡點的移動能力進行定義，其劃定鄰域半徑R作為軌跡點移動能力計算范圍，但對于存在乒乓切換和信號漂移的手機信令而言，鄰域半徑R將軌跡的時間序列切割為多個軌跡段，而移動能力最小的軌跡段不能有效反映該連續時間序列的軌跡段移動能力。

本文網格軌跡N經過時空聯結后形成網格簇序列{Traj1，Traj2，…，Trajf}（f為聯結后的網格簇總數）。每個網格簇的移動能力通過自身軌跡與該網格簇和前后網格簇之間轉移軌跡的影響，以及網格簇內的軌跡來定義。假設網格簇自身軌跡為Traji={pm，pm+1，…，pn}，該網格簇與前后網格簇之間轉移軌跡為dis(pm－1，pm)與dis(pn，pn+1)之和，轉移時間為travel time(pm－1，pm)與travel time(pn，pn+1)之和，網格簇總時間為移動時間與停留時間之和。

定義3（網格簇移動時間）網格簇移動時間是當前網格簇Traji={pm，pm+1，…，pn}與前后網格簇轉移時間以及網格簇內單網格與單網格、單網格與停留網格的轉移時間之和，計算公式如下：

其中，v為當前網格簇內單網格與停留網格的總個數。

定義4（網格簇總時間）網格簇總時間是從網格簇Traji={pm，pm+1，…，pn}的前一點pm-1到網格簇的后一點pn+1所經歷的時間，計算公式如下：

定義5（網格移動能力）網格移動能力是前后網格簇和當前網格簇的轉移軌跡與當前網格簇內軌跡自身曲折性和移動時間比例的乘積，計算公式如下：

由式（6）可知：軌跡起始點與目的地之間軌跡越曲折，則網格移動能力MA 越小；整個軌跡運動中移動時間占總時間比例越小，則網格移動能力MA 也越小，與經驗認知相符。

根據數據域理論，空間中每個數據點都會受其周圍數據點的影響，且該影響隨距離的增加而減弱。由于軌跡停留區域內數據點聚集程度明顯高于移動區域內數據點，因此停留區域內每個點所受影響也高于移動區域內數據點［22］。與以軌跡點為研究單元，并以軌跡點與其周圍軌跡點間的影響定義軌跡點密度不同，本文以網格簇軌跡為研究單元，結合數據域理論，通過時間序列前后網格簇及網格簇自身影響定義網格簇時空密度函數MAST，計算公式如下：

其中，σMA和σdis分別為移動能力權重函數和距離權重函數的標準偏差，ST(pk)=tk+1-tk為網格簇轉移軌跡與簇內軌跡網格間的停留時間。

由式（7）可知：網格軌跡移動能力MA 越小，則網格間距離越小，停留時間越長；網格簇時空密度函數MAST 越大，則其為停留簇的可能性越大。與以往時空密度的定義不同［21-22］，本文以網格簇為對象，利用時空聯結形成網格簇來實現簇內軌跡的時空連貫性，同時簇與簇之間的轉移距離與轉移時間進一步保證手機信令的時間序列性不被空間鄰域半徑所割裂，并與軌跡曲折性相結合，充分應用手機信令特性，增強了時空密度判斷的合理性。

通常采用速度來判斷是否為移動和停留狀態，但是由于手機信令的固有特性，用戶定位基站并不代表用戶的實際位置，且基站覆蓋范圍較廣，存在覆蓋區域相互重疊的情況，因此基站內停留時間較實際時間偏大，基站間轉移時間較實際時間偏小，從而造成基站間移動速度較實際速度偏大。同時，因為存在噪聲，乒乓切換和信息漂移的速度均較大，所以速度波動幅度較大。將各網格簇的平均移動速度與其時空密度函數值進行對比，結果如圖3 所示。

圖3 網格簇的平均移動速度和時空密度Fig.3 Average moving speed and space-time density of grid clusters

由圖3 可以看出：移動速度大的區域與許多速度小的區域交織，說明移動過程中存在個別減速或擁堵時段；網格簇的移動速度與時空密度的變化相反，移動速度大的網格簇時空密度小，而移動速度小的網格簇時空密度大，說明時空密度能夠反映移動和停留的特性；時空密度較移動速度更穩定，對移動和停留更有區分性，因此其更適合作為判斷移動和停留狀態的指標。

時空密度函數參數包括密度函數中的σMA、σdis以及時空密度函數MAST 的最小密度閾值。對于σMA和σdis的選取，本文通過自行開發的基站采集APP 采集真實手機信令軌跡數據來驗證參數選取的可靠性。將不同的σMA和σdis進行組合計算得到各網格簇的時空密度函數MAST，通過其數據概率密度分布自動確定最小密度閾值，篩選得到大于時空密度函數MAST 的最小密度閾值的網格簇個數，對比不同σMA和σdis組合下停留簇的變化情況，結果如圖4 所示。

圖4 不同σMA和σdis組合下網格簇時空密度的概率密度分布Fig.4 Probability density distribution of space-time density of grid clusters with different combinations of σMAand σdis

由圖4 可以看出，σMA的取值為0.1～0.9，σdis的取值為0.5～50，不同σMA和σdis組合下網格簇的MAST 的概率密度分布基本相同。由此可知，σMA和σdis的取值對MAST 的概率密度分布影響較小，且σMA對MAST 的影響大于σdis。此外，MAST 概率密度分布曲線有兩個峰值，較大峰值為用戶職住地的時空密度，較小峰值包含移動簇和短時間的停留簇。根據經驗可知，用戶在工作地和居住地的活動時間通常占一天內大部分時間，其時空密度遠大于移動點與短時間停留區域，因此兩個峰值之間存在一段MAST=0 的分布。因此，可先采用K-means 算法將MAST 聚類為職住地、常規移動點和停留點兩大類，再將MAST 密度小的一類（常規移動點和停留點）用K-means 算法再聚類為兩類，將這兩類的分界值作為MAST 最小密度閾值以區分移動簇和停留簇。綜上所述，將MAST 用K-means 算法聚類為A、B、C和D 4類，不同σMA和σdis組合下MAST的聚類效果如圖5 所示（彩色效果參見《計算機工程》官網HTML 版）。

圖5 不同σMA和σdis組合下網格簇時空密度的K-means 聚類分布Fig.5 K-means clustering distribution of grid cluster space-time density under different combinations of σMAand σdis

由圖5 可以看出，在不同σMA和σdis組合下網格簇MAST 的K-means 聚類分布中，移動簇和停留簇的分布基本相同，均識別出7 個停留區域，說明參數σMA和σdis的選取對停留區域的識別影響較小，通過K-means 聚類算法自動選取MAST 閾值得到的識別結果具有較好的穩定性。

2.4 網格簇時空聚類

從網格簇中識別出移動簇與停留簇后，將時空距離較近的簇進行時空聚類，以錨固停留區域并避免停留區域被噪聲分割。

根據停留和出行的定義，一次停留時間需大于5 min，一次出行時間需大于5 min 且出行距離大于500 m。本文定義一次出行起訖點簇最小覆蓋圓的圓心距離需大于500 m 且兩圓相離，設計時空聚類流程如下：

1）計算各停留簇最小覆蓋圓的圓心c、半徑r以及簇停留時間st=tn－tm，即Traji={pm，pm+1，…，pn}中軌跡點首尾時間戳的差值。

2）若當前停留簇Traji={pm，pm+1，…，pn}與下一個停留簇Traji+1={px，px+1，…，py}首尾時間間隔tx-tn>5 min，當前停留簇Traji最小覆蓋圓與下一個停留簇Traji+1最小覆蓋圓相離，且其圓心ci與ci+1的距離大于500 m，滿足出行條件，則轉到步驟3；否則轉到步驟5。

3）若當前停留簇Traji的停留時間st(i)>5 min，滿足停留條件，則將停留簇Traji標記為停留區域，下一個停留簇成為當前停留簇；否則轉到步驟4。

4）若當前停留簇Traji的停留時間不滿足停留條件，則將該停留簇標記為延誤區域。

5）將下一個停留簇Traji+1加入當前停留簇Traji，合成Traji={pm，pm+1，…，pn，px，px+1，…，py}，重新計算最小覆蓋圓的圓心ci、半徑ri以及簇停留時間st(i)，轉到步驟2。

本文提出的網格簇時空聚類算法，考慮了停留簇之間時間間隔與空間的關系，是一種漸進的聚類算法［23］。該算法將一定區域的小空間簇進行聚合有利于錨固有意義的停留區域。

3 算例與結果分析

本文采用自行開發的基站采集APP 采集志愿者的軌跡數據（包括時間戳、基站經緯度、移動或停留標簽數據）以驗證算法的識別效果。基站采集頻率設置為5 s 采集一次（由于各品牌手機硬件不同，因此，實際采集頻率大于該值）。記錄志愿者一天24 h 的軌跡數據，得到10 582 條數據，其中包括5 個停留區域與3 次交通換乘記錄，由于換乘時間均大于5 min，因此3 個換乘站點在軌跡識別結果中均視為停留區域。將10 582 條記錄進行網格化，并與同網格數據合并，只保留首尾兩條數據，將數據簡化為846 條，并通過時空聯結操作將846 條網格數據劃分為時空緊密相連的144 簇。表1 為時空聯結生成的網格簇軌跡記錄。以簇2 為例，將序號4 和序號5 記為組1，序號7 和序號8 記為組2，序號6 記為組3，序號9 和序號10 記為組4，組1 和組2、組3 和組4 分別表示簇2 內兩個地理位置相同和兩個地理位置相鄰且時間間隔小于15 min 的網格簇軌跡記錄，上述各組區域相互交織，共同形成簇2。

表1 時空聯結生成的網格簇軌跡記錄Table 1 Track record of grid clusters generated by space-time connections

從空間識別有效性和時間識別有效性將本文算法與改進DBSCAN 時空密度算法［21］（以下稱為文獻［21］算法）得到的停留點識別結果進行對比，采用停留點個數、停留時間的準確率P、召回率R以及準確率和召回率的加權平均值F1 來評價本文方法的有效性，并通過對比兩種算法的耗時以評價算法效率。文獻［21］算法參數設置為：σMA=0.5，σdis=0.5，Nap=51。本文算法參數設置為：σMA=0.5，σdis=0.5。兩種算法得到的停留點識別結果如表2 所示。可以看出：文獻［21］算法檢測到的停留區域個數在準確率和召回率上均明顯低于本文算法，本文算法在空間上識別有效性更高；在時間的識別有效性方面，文獻［21］算法檢測到的停留時間等于與實際相符的停留時間且遠小于標記停留區域的停留時間，停留區域個數多于標記的停留區域個數，說明文獻［21］算法難以檢測到停留區域的邊緣點，且在噪聲影響下將一個大的時空停留區域切分成幾個小區域。

表2 兩種算法的停留點識別結果Table 2 Recognition results of stop points of the two algorithms

在時間效率方面，文獻［21］算法的時空密度計算次數為Nap×N（N為軌跡點的個數），本文算法的計算次數為網格軌跡時空聯結后所得簇數S，文獻［21］算法作為DBSCAN 聚類算法的變型，其復雜度為O(N2)。在計算效率方面，本文算法的耗時遠低于文獻［21］算法，計算效率更高，可直接應用于大樣本的手機信令數據計算，而文獻［21］算法無法應用于大樣本數據的計算。

圖6 為標記的停留區域與不同算法檢測到的停留區域（彩色效果參見《計算機工程》官網HTML版）。其中，標記的停留區域stay4 被文獻［21］算法分割為stay7、stay8（被stay9 覆蓋）、stay9、stay10 和stay11 共5 個小區域，而本文算法檢測到的停留區域與標記的停留區域stay4 地理位置一致，說明本文算法能將時空緊密的軌跡點聯結為一個整體，并有效識別出停留區域，在較大程度上消除噪聲以及采樣不均造成的空間不確定性對停留點識別的影響，從而錨固住停留區域，不易使停留時空區域被分割。此外本文算法中時空移動能力相較文獻［21］算法具有更好的區分性與穩定性，可使網格簇的時空密度測定更準確。

圖6 標記的停留區域與不同算法檢測的停留區域Fig.6 Marked residence area and residence area detected by different algorithms

由上述分析可知，文獻［21］算法的可靠性除了選定σMA和σdis外，還與采樣間隔及Nap 參數選取相關，而參數選取對本文算法影響較小，可根據時空密度分布自動選取閾值，因此本文算法的穩定性和可靠性較文獻［21］算法更優。結合圖6 中數據質量來看，由于本文實驗采集的出行軌跡較密集，而現實中手機信令頻率遠小于實驗采集頻率，合并連續同位置的網格與現實中信令數據的觸發更接近，且本文基于網格簇空間密度劃分軌跡點，較基于軌跡點密度劃分軌跡點在空間上的可行性更強，因此本文算法更適用于現實中手機信令軌跡點識別。

4 結束語

手機信令因具有數據采樣率不均、定位精度低與空間不確定性，給網絡軌跡點分析、出行鏈提取等基礎性研究帶來困難，而現有消除空間不確定性的預處理算法缺乏從時空角度綜合度量軌跡點移動與停留情況的指標。本文提出一種時空密度軌跡點識別算法，根據出行定義和振蕩噪聲特征通過時空聯結識別潛在停留簇，以消除空間不確定性對后續聚類的影響，同時結合軌跡的曲折性、移動時間和停留時間重新定義網格簇內軌跡點的移動能力，基于前后網格簇間與簇內軌跡點間的時空關系和軌跡簇移動能力構建時空密度指標進行軌跡點停留區域分析。實驗結果表明，與改進DBSCAN 算法相比，該算法識別精度和時效性更高。由于軌跡點識別結果是具有停留屬性的軌跡點集群，其中包括真實活動區域、候車區域以及交通擁堵區域等，因此后續將對停留區域的真實活動屬性進行研究，以挖掘個體活動規律，進一步提高識別精度和效率。