一種聚類隱馬爾可夫模型的時空軌跡預測算法

孫 紅，陳 鎖

1(上海理工大學，上海 200093) 2(上?，F代光學系統重點實驗室，上海 200093)

1 引 言

近年來，由于大數據時代的來臨，以及物聯網技術的廣泛應用，各種移動終端的大量普及，大量GPS數據隨之產生，各種軌跡數據也越來越容易獲得，對GPS軌跡數據的研究也越來越多.大部分終端設備都已裝有GPS模塊，這些終端設備因而可以記錄GPS位置序列，這些序列包含經緯度，時間戳等信息.時空軌跡(Trajectory)是移動對象的位置和時間的記錄序列.作為一種重要的時空對象數據類型和信息源，時空軌跡的應用范圍涵蓋了人類行為、交通物流、應急疏散管理、動物習性和市場營銷等諸多方面.時空數據的急劇增加，加之時空數據處理更為復雜，使數據處理任務日趨繁重的形勢更加嚴峻.因此，尋找有效的時空數據處理方法具有十分重要的意義.

時空大數據與非空間數據相比，具有空間性、時間性、多維性、海量性、復雜性等特點.關于時空大數據，目前存在一些技術前沿領域，如云計算方法和挖掘技術.而預測移動對象的運動趨勢有著廣泛的應用.比如能夠有效地在軌跡數據中發現隱藏的有價值的信息，如頻繁路徑、個性化興趣點以及移動意圖等.這些結果將對未來的商業應用以及管理活動提供有力的信息保障.

本文通過分析時空軌跡數據，提出一種根據移動對象經過的軌跡，預測移動對象下一個軌跡點，以及未來的終點位置的模型.

2 研究背景

對于移動對象軌跡位置預測，近年來主要產生了兩種類型.一種是對歷史軌跡頻繁模式進行挖掘進而進行預測.另一種是基于馬爾可夫(Markov)模型，對歷史軌跡進行建模并預測.

軌跡頻繁模式挖掘算法，目前主要分為兩種：Apriori算法和PrefixSpan算法.一些研究將它們應用到軌跡位置預測的建模中.一種基于廣度優先搜索的Apriori算法結合自定義的運動函數和移動軌跡模式，堪稱一種創新型的方法，可以參考Jeung等人[10]的相關文章.還有一種改進版的Apriori算法，可參考Morzy[11]的相關文章.另一種主流算法是PrefixSpan算法，文獻[12]提出一種改進型的PrefixSpan用來發現運動軌跡的歷史頻繁項并生成預測模型.此外，Monreale等人[2]提出的WhereNext方法也不失為一種好的方法，該算法從歷史軌跡中提取具有幾種不同運動行為的運動模式，進而生成預測模型.

馬爾可夫模型對軌跡序列建模中，有一些比較新的成果，比如Figueiredo等人[7]，他們使用不穩定，瞬態且時間多樣化的軌跡序列，構造出一種非參數模型用來對軌跡位置等做出預測.文獻[9]使用了一種HITS-based模型來挖掘用戶軌跡興趣點，側重于分析軌跡全局的運動模式.

總的來說，現有工作大多是依據大量移動對象軌跡的重合部分，挖掘出軌跡相似行為，再用來做預測.這種方式的缺點是：著重考慮表面發生的行為，對于事物內在的隱含信息的考慮有所欠缺.

由于移動物體產生的軌跡數據量比較大，而且這樣的數據還在不斷地增長，因而對于以上兩種方法，模型的訓練成本非常大.本文基于此原因，提出一種基于聚類分區域的隱馬爾可夫模型，經過大量GPS數據的訓練，得出最佳模型參數，最后通過隱狀態之間的轉移矩陣預測移動對象軌跡最終位置.

3 聚類分區域的隱馬爾可夫模型

3.1 隱馬爾可夫模型(HMM)

馬爾可夫模型中，狀態是可見的，因此其又被稱為可視馬爾可夫模型(visible Markov model，VMM)，狀態可見，在一定程度上限制了模型的適應性.而在隱馬爾可夫模型(HMM)中，每個狀態是不可見的，可見的是各個狀態產生的事件的概率分布，狀態之間的轉換也是不可觀察的，可觀察的事件是由狀態產生的，每個狀態產生一個特殊分布的事件集合.圖1描述了隱馬爾可夫模型的基本原理.

圖1 隱馬爾可夫模型Fig.1 Hidden Markov model structure

一個HMM由如下幾個部分組成：

1)模型中隱狀態的數目N；

2)隱狀態可能輸出的不同可觀察符號的數目M

3)狀態轉移矩陣A={aij}，其中，

aij=P(qt=sj|qt-1=si)，1≤i，j≤N

aij≥0

4)可觀測事件概率分布矩陣B={bj(k)}，其中，

bj(k)=P(Ot=vk|qt=sj)，1≤j≤N；1≤k≤M

bj(k)≥0

5)初始狀態概率分布π={πi}，其中，

πi=P(q1=si)，1≤i≤N

πi≥0

一個HMM一般來說可以用以一個三元組來描述，μ=(A，B，π)，A為狀態轉移概率，B為符號發射概率，π為初始狀態概率分布.當一些現象產生是由于某些事物本質變化做決定時，用HMM進行建模一般是非常有用的.因為HMM能夠“深刻理解”事物的變化規律.HMM三個基本問題及其解決方法如下：

1)給定一個模型μ=(A，B，π)，如何快速計算出觀察序列O=O1O2…OT的概率，即P(O|μ).

解決方法：前向算法、后向算法.

2)給定一個模型μ=(A，B，π)，和一個已經產生的觀察序列O=O1O2…OT，如何找出最可能產生觀察序列O的狀態序列Q=q1q2…qT.

解決方法：維特比算法.

3)給出大量觀察序列集，如何構造出一個合適的隱馬爾可夫模型μ=(A，B，π).解決方法：Baum_Welch算法

HMM在自然語言處理研究中有著非常廣泛的應用，尤其是在語音識別領域，一段語音可以作為觀測序列，語音對應的文字作為隱狀態，語音識別的首要任務就是求出最有可能發出該語音序列的隱狀態，將語音恢復成文字.通過一部分可觀測序列，可以預測對應的隱狀態序列，基于此，本文將HMM模型移植到時空軌跡序列預測的問題中.將大樣本時空軌跡序列，作為HMM的可觀測序列，將局部的時空坐標點看成是附近某個隱狀態所發射出的觀測值，該隱狀態可以視為該局部區域的地標點.如圖2所示.

圖2 隱狀態及觀測值Fig.2 Hidden state and visible state

圖中黑色三角形即為隱狀態，表示城市內的多個地標，這些地標可以看成是人群最愛去的地點的集合.以A點為例，周圍有1，2，3，4，5，6幾個GPS坐標點，它們分別是四條GPS軌跡中的點，因為都處于A點附近，因此可認為它們是A狀態發射出的觀測點，即它們服從A的發射概率分布.本文的任務之一就是要建立這樣的HMM模型，模型訓練完成之后，輸入要預測終點的一段時空軌跡序列，即為觀測序列，然后通過相應的算法求解出最佳隱狀態序列，通過隱狀態序列最后一個狀態以及轉移矩陣，計算出概率最大的下一個隱狀態，再根據該隱狀態的發射概率分布計算出最有可能的發射值即為預測的終點.整個流程如圖3所示.

3.2 聚類分區域的隱馬爾可夫模型——CSHMM

HMM模型常用來進行語音識別，一段語音的音節數量是很少的，一般是幾十個，這樣的數量級，采用維特比算法來求解隱狀態，即推斷語音對應的文字或單詞，算法的執行時間是非?？焖俚?其原因是維特比算法的時間復雜度是O(N2T)，其中N為發射出來的可觀測點的數目，T為狀態序列的長度.由于車輛一般都是往一個大致的方向行駛(一般不會存在繞著整個城市或局部區域反復行駛)，因而一條完整的行駛路徑所穿過的地標區域個數是比較少的，因而用維特比算法也是很快可以計算出最佳隱狀態序列.然而整個城市的所有地標區域是非常多的，比如商場，運動場，小公園，小吃街等，實際上任何地方都可能成為車輛停止的終點.因為時空序列終點的預測精度只要精確到一定的范圍內即可，因此所有隱狀態的數量就非常多了，可能會有幾萬甚至幾十萬個.如果針對整個城市建立一個含有如此龐大數量隱狀態的隱馬爾可夫模型(尤其是北上廣深這樣的特大城市)，那么在用維特比算法求解隱狀態序列時，遞歸算法將會非常耗時，甚至很可能導致程序直接奔潰.本文針對該情形，設計了一種聚類分區域的隱馬爾可夫模型(Clustering Subregion Hidden Markov Model，即CSHMM).

圖3 算法流程圖Fig.3 Algorithm flow chart

3.2.1 產生隱狀態集

建模的第一步，是要確定所有的隱狀態集合，在本文中即為確定整個區域中的地標點集合.

一般來說GPS經緯度坐標數值精確到小數點后4位即可明顯的區分兩個坐標點為兩個不同的地點(相差幾十米至150米的兩個地點本文當作同一地點)，統計訓練數據集中的所有不同地點出現的頻率fi.設置一個閾值r，將fi>r的所有地點標記為地標.每個地標使用其頻率作為權重.

其次，將得到的所有地標點集，根據其權重進行聚類，聚類的原則是：權重越高的地段所包含的區域面積越小，以使多個子區域包含的數量級相當.在此我們采用坐標點的數量來代替區域面積，可以對此建立一個線性模型(也可以是非線性模型)來確定每個區域內的大致坐標點數目.本文采用一種改進型的kmean——kmeans_cs如下函數模型來判定每個地標周圍的坐標點是否屬于該地標的周邊區域：

設o為一個聚類中心點，p為需要歸類的點，n為o類樣本集中包含p點的所有軌跡數量，o和p之間的距離D(o，p)采用如下定義：

D(o，q)=e‖o，p‖·‖o，p‖

(1)

其中||o，p||為o，p兩點的歐氏距離.由公式(1)可以看出兩個坐標點的距離越近，則會有一個大于1的權重乘進去使其增大.如果距離越大，那么乘上的權重會以指數增大，使得距離增加的更快.這樣會使得聚類的結果中面積越小的子區域包含的點越多，面積大的區域包含的點越少.改進后的聚類算法可參考下文中2.2.2節中的聚類算法.

圖4大致描繪出映射建立結束后的圖形.

圖4 分區域后的隱狀態Fig.4 Hidden state after subregion

圖5中的黑色三角形即為隱狀態，周圍的點即為隱狀態的發射符號集.

直觀表示如圖5所示.

圖5 隱狀態和觀測值的映射Fig.5 Mapping of hidden states and observation values

如圖4所示，首先通過聚類算法將城市中的所有地標點進行聚類，以此將全局區域分成若干個小區域，然后針對每個小區域，進行HMM建模.模型建立完成以后，根據輸入的待預測的部分GPS軌跡，首先進行區域匹配，選出最佳HMM模型，然后使用維特比算法針對該模型進行預測.

3.2.2 聚類劃分子區域

由于地標點數目太多，因而我們需要劃分子區域，使得每個子區域內的地標點根據其權重均勻分布.劃分的原則是：子區域之間GPS軌跡跨越盡可能少.本文針對GPS軌跡的變化規律設計了一種改進的k-means聚類算法，傳統的k-means算法采用歐氏距離來判斷點和點之間的類似程度，此方式會打破GPS軌跡坐標的連貫性，即會大量出現同一條軌跡存在于多個子區域中，這是不符合要求的，因為本文需要在每個子區域單獨求解HMM模型，每個模型模擬了一個特定子區域GPS軌跡模式.基于該目的，本文提出一種能夠按照GPS軌跡來聚類的算法，該算法改進了點和點之間的距離公式：

(2)

設o為一個聚類中心點，p為需要歸類的點，n為o類樣本集中包含p點的所有軌跡數量，o和p之間的距離D(o，p)采用如下定義：

D(o，p)=tanh(n)×‖o，p‖

(3)

其中‖o，p‖為o，p兩點的歐氏距離.采用tanh函數來作為權重函數.不采用sigmoid是因為sigmoid自變量取值為0時函數值為0.5，而這里想要的是在區域p內含有q點的軌跡數目為0時，權重就變為0.

改進后的聚類算法如下：

1)任選K個初始聚類中心Z1(1)，Z2(1)，Z3(1)，…Zk(1)；

根據樣本間的相互距離，將樣本分配到K個聚類區域，即若Min{d(Zi(k)，X)，i=1，2，…K} = d(Zj(k)，X) = Dj(k)，則X∈Sj(k)，其中，k為迭代序號.

2)更新各個聚類中心Zj(k+1)，其中j=1，2，3...K；

3)若Zj(k+1)≠Zj(k)，則跳轉到步驟(2)，將各個樣本重新分類，如此迭代；如果Zj(k+1)=Zj(k)，j=1，2，…K，算法收斂，計算完畢.

3.3 分區域進行HMM的模型訓練

對于每個子區域，其所有的地標數據已經確定，每個地標周圍的稱為發射符號集的GPS坐標點也已經確定.已知的信息量比較充足，比如可以根據每個GPS坐標的出現頻率，計算出發射概率矩陣B0.也可以根據大量軌跡來計算地標與地標之間的轉移概率，從而計算出轉移矩陣A0.

Baum-Welch算法是一種EM算法，專門用來訓練HMM模型.EM過程保證算法一定能夠收斂到一個局部最優點，但是它不能保證找到全局的最優點.所以Baum-Welch算法不能保證找到全局的最優模型.但是由于我們已經知道了大致的解：狀態轉移矩陣A0和發射概率矩陣B0，如果將它們作為迭代的初始值，那么會有很大的概率收斂到全局最優點，并且訓練的時間也大大減小.本文采用的Baum-Welch算法.

對于給定的隱馬爾可夫模型(HMM)μ和已知的觀察序列O=O1O2…OT，假設在t時刻該模型位于狀態si，在時間t+1時刻模型位于狀態sj，此概率ξt(i，j)可由如下公式計算獲得：

(4)

給定HMM的參數μ和觀察序列O=O1O2…OT，在時間t位于狀態si的概率γt(i)有下面的公式計算獲得：

(5)

而μ的參數可以有下面的公式重新估計：

(6)

(7)

(8)

步驟1. 隨機給參數πi，aij，bj(k)賦值，并使其滿足如下約束：

由此得到新的模型μ0，令i=0，并用期望最大化算法(EM)重新估計模型的各個參數；

步驟2. EM計算：

E步：由模型μi根據公式(4)和公式(5)計算出ξt(i，j)和γt(i)；

M步：用E步得到的期望值，根據公式(6)，公式(7)，公式(8)重新估計參數πi，aij，bj(k)，得到模型μi+1；

步驟3. 循環計算：

令i=i+1.重復執行EM計算，直到πi，aij，bj(k)收斂.

3.4 使用模型進行預測

對于給定的部分GPS軌跡序列，根據最后一個坐標的位置，確定相應的子區域，選取相應子區域的HMM模型.將不在該子區域內的坐標點切除.使用靠近尾部的坐標點作為可觀測序列.問題就變成了HMM模型的第二個問題.因此本文采用維特比算法來求隱狀態序列.算法如下：

步驟1. 初始化：

δ1(i)=πibi(O1)，1≤i≤Nφ1(i)=0

步驟2. 歸納計算：

記憶回退路徑：

步驟3. 終結：

步驟4. 路徑回溯：

求出隱狀態序列后，根據轉移矩陣A=[aij]和最后一個隱狀態qlast，求出下一個概率最大的可觀測坐標點Pnext.求法如下：

(9)

其中last為當前時刻，next為下一時刻.

4 實驗與分析

本文使用的數據集為T-Drive Trajectory data sample，來自微軟亞洲研究院(MSR)，一共有10357個北京出租車在一周內的GPS軌跡，每條軌跡提取出包含經緯度和時間戳信息的子信息.實驗選區75%的數據作為訓練數據集，25%的數據作為測試數據集.代碼使用python語言，硬件配置為 Intel(R) Core(TM) i5-3320 2.60GHZ，12.0GB內存，320GB硬盤的計算機.

實驗結果采用準確率，召回率為評估標準.

4.1 本文模型使用kmeans和改進型kmeans結果對比

結果如圖6所示(準確率和召回率為改進kmeans后的結果，準確率2和召回率2為不使用改進的kmeanss).

圖6 kmeans和改進型kmeans對模型效果的影響Fig.6 Effect of kmeans and improved kmeans on the model effect

結果分析：改進kmeans后的算法導致模型準確率和召回率提高了15%左右，說明改進后的kmeans認為同一個軌跡中的點更應該屬于同一個聚類，因為這些點之間是有關聯的，有關聯的點之間的跳轉，直接導致轉移概率的增大，因而會有更優秀的準確率.

4.2 本文模型和其他模型對比

結果如圖7所示(準確率和召回率為改進kmeans后的結果，準確率3和召回率3為使用Apriori算法的結果，橫坐標單位：萬條數據，縱坐標單位：百分比).

圖7 本文模型和其他模型性能對比Fig.7 Comparison of the performance of this model and other models

結果分析：改進kmeans后的算法導致模型準確率和召回率提高了8%左右.對比傳統關聯模式挖掘的算法，本文算法具有一定的優勢.說明基于轉移概率的方式與挖掘時空序列軌跡相似度的方式相比，前者將一些影響軌跡趨勢的別的特征也考慮進了模型中(即使沒有提取相關特征)，因而會提高準確度.

5 結束語

本文基于隱馬爾可夫模型，提出了一種預測移動對象下一個位置的方法—基于聚類分區域的時空軌跡預測模型，該方法用以解決時空軌跡序列的預測問題.該模型首先通過聚類將一片區域內的時空序列分成多個小區域，每個小區域內再通過聚類確定多個隱狀態和發射序列，然后針對每個小區域進行隱馬爾可夫模型的訓練得出最終模型.預測時通過已知的時空序列，找到對應的區域模型，通過維特比算法計算出最佳隱狀態序列，再結合轉移矩陣做出下一個軌跡點的預測，實驗表明該方法在一定程度上提高了預測準確率.

未來的工作中，可以使用機器學習或者深度學習的方式進一步提高區域的劃分算法和軌跡聚類算法.