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

陳 略，熊 宸，蔡 銘

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

0 概述

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

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

1 相關工作

1.1 手機信令的特性

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

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

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

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

1.2 手機信令預處理

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

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

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

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

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

1.3 停留點識別

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

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

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

基于軌跡形狀［11］用定義的最小覆蓋圓沿時間軸將軌跡點劃分到簇中，當軌跡點超過最小覆蓋圓時判斷簇的停留時間，若大于停留時間閾值則判斷為停留簇，若小于停留時間閾值則排除簇中時間軸最靠前的軌跡點，并從簇中第二個軌跡點開始循環(huán)?；谛螤畹木垲愃惴ú皇芫垲愴樞虻挠绊?，較基于時間聚類的算法能更有效地劃分距離近的簇，但容易受到噪聲影響而分為多個小簇，且該算法受停留時間閾值的影響較大。

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

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

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

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

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

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

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

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

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

2 本文算法

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

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

2.1 軌跡點網(wǎng)格化

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

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

2.2 網(wǎng)格軌跡時空聯(lián)結

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

算法1網(wǎng)格軌跡時空聯(lián)結算法

上述算法具體步驟如下：

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

2）按時間序列輸入當前軌跡點記錄(i，xi，yi，ti)，若下一條記錄(i+1，xi+1，yi+1，ti+1)為相鄰網(wǎng)格，則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中，當前記錄變?yōu)橄乱粭l記錄(i+1，xi+1，yi+1，ti+1)，即i=i+1，然后再從步驟2 開始；否則轉(zhuǎn)到步驟4。若下一條記錄所在網(wǎng)格與當前網(wǎng)格不相鄰，則轉(zhuǎn)到步驟3。

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

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

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

算法2刪除與循環(huán)聯(lián)結算法

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

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

2.3 網(wǎng)格時空密度

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

圖2為兩種軌跡的移動情形。可以看出在同一個起始點和目的地的軌跡段中，移動的軌跡偏向進行直線運動（見軌跡1），而停留或徘徊的軌跡更偏向進行曲線運動（見軌跡2）。文獻［21］根據(jù)圖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 的定義僅考慮整個軌跡段的曲折性，而長時間在職住地或者基站覆蓋范圍較大的區(qū)域停留時，與其他基站的切換較少但停留時間較長，如果只考慮曲折性，那么在停留區(qū)域位置不變的點的移動能力為1，其移動能力反而會異常大。因此，本文綜合考慮軌跡的曲折性、移動時間和停留時間后，對網(wǎng)格時空移動能力進行重新定義。文獻［21］針對核心點通過該點所在軌跡段的MA 及其周圍網(wǎng)絡接入點（Nap）對其產(chǎn)生的影響來定義，但其利用周圍軌跡點個數(shù)劃定移動能力的計算范圍不適用于采樣頻率不均的手機信令數(shù)據(jù)。文獻［22］基于軌跡點鄰域半徑R內(nèi)包含所有軌跡段的最小值對軌跡點的移動能力進行定義，其劃定鄰域半徑R作為軌跡點移動能力計算范圍，但對于存在乒乓切換和信號漂移的手機信令而言，鄰域半徑R將軌跡的時間序列切割為多個軌跡段，而移動能力最小的軌跡段不能有效反映該連續(xù)時間序列的軌跡段移動能力。

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

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

其中，v為當前網(wǎng)格簇內(nèi)單網(wǎng)格與停留網(wǎng)格的總個數(shù)。

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

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

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

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

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

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

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

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

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

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

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

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

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

2.4 網(wǎng)格簇時空聚類

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

根據(jù)停留和出行的定義，一次停留時間需大于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，滿足出行條件，則轉(zhuǎn)到步驟3；否則轉(zhuǎn)到步驟5。

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

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

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

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

3 算例與結果分析

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

表1 時空聯(lián)結生成的網(wǎng)格簇軌跡記錄Table 1 Track record of grid clusters generated by space-time connections

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

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

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

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

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

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

4 結束語

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