基于主成分回歸算法的城市客流聚集風(fēng)險(xiǎn)預(yù)測(cè)

王聚全， 王 偉， 馬慧民， 楊 博， 杜 渂

1.迪愛(ài)斯信息技術(shù)股份有限公司，上海200032

2.電信科學(xué)技術(shù)第一研究所有限公司，上海200032

3.上海市公安局科技處，上海200042

4.上海市北高新股份有限公司，上海200436

伴隨中國(guó)經(jīng)濟(jì)的騰飛，大型活動(dòng)逐漸成為促進(jìn)文化交流和經(jīng)濟(jì)繁榮發(fā)展的重要載體.諸如上海、北京等規(guī)模較大城市，各種大型活動(dòng)（如節(jié)假日歡慶、體育賽事等）的規(guī)模和數(shù)量與日俱增，隨之引發(fā)的各種群體性安全問(wèn)題層出不窮.據(jù)中國(guó)社科院公布的《社會(huì)藍(lán)皮書(shū)》數(shù)據(jù)顯示，中國(guó)發(fā)生的群體性事件從1993年至2005年約增加10 倍，且近幾年依然呈現(xiàn)上升趨勢(shì)[1-2]，尤其是2014年元旦跨年夜發(fā)生在上海外灘的踩踏事件再一次給政府應(yīng)急管理敲響了警鐘.

對(duì)以往發(fā)生的群體性事件進(jìn)行分析后不難發(fā)現(xiàn)，客流密度過(guò)大、過(guò)度擁擠是事故發(fā)生的根本原因[3-4].在大城市，移動(dòng)網(wǎng)絡(luò)基本實(shí)現(xiàn)了全覆蓋，且手機(jī)的人群使用量也達(dá)到了較高水平[5].因此，如何有效利用現(xiàn)有數(shù)據(jù)預(yù)測(cè)客流量的變化情況已逐漸成為業(yè)界的重點(diǎn)研究方向，目前客流預(yù)測(cè)方法大致可分為兩類(lèi)：

一類(lèi)是基于時(shí)間序列的模型預(yù)測(cè)法[6-7]，代表算法有差分自回歸移動(dòng)平均（autoregressive integrated moving average, ARIMA）模型.該類(lèi)算法具有模型簡(jiǎn)單、不需要借助輸入變量的優(yōu)點(diǎn)，但不能捕獲非線(xiàn)性關(guān)系且有輸入數(shù)據(jù)平穩(wěn)或差分后平穩(wěn)的要求.

另一類(lèi)是基于降維的模型預(yù)測(cè)法[8]，代表算法有主成分分析（principal component analysis, PCA）法.該類(lèi)算法的優(yōu)點(diǎn)是在損失較少的情況下能從眾多信息中抽取關(guān)鍵信息，缺點(diǎn)是需要保證主成分因子符號(hào)為非負(fù).

綜合以上兩類(lèi)算法的優(yōu)缺點(diǎn)以及通常利用運(yùn)營(yíng)商提供的手機(jī)用戶(hù)數(shù)據(jù)這兩方面因素，本文采用主成分分析、回歸分析和最小二乘法相結(jié)合的數(shù)學(xué)模型對(duì)特定區(qū)域的客流聚集人數(shù)和趨勢(shì)進(jìn)行預(yù)判，與ARIMA 算法相比在預(yù)測(cè)準(zhǔn)確度方面有較大提升.

1 算法模型理論

1.1 ARIMA 模型

時(shí)間序列的分析模型包括以下幾種：自回歸[9]（auto-regressive,AR）、移動(dòng)平均[10]（moving average, MA）、自回歸移動(dòng)平均[11]（auto-regressive moving average, ARMA）、ARIMA 等模型，其中ARIMA 模型作為近年來(lái)的最新研究成果，受到了工商界和學(xué)術(shù)界的廣泛關(guān)注，其核心公式為

式中，?[.]表示AR 的系數(shù)，e[.]表示MA的系數(shù).ARIMA 模型要求輸入的時(shí)間序列是平穩(wěn)時(shí)間序列或者經(jīng)過(guò)差分后的平穩(wěn)時(shí)間序列，檢驗(yàn)時(shí)間序列是否平穩(wěn)的技術(shù)手段是使用時(shí)間序列的自相關(guān)系數(shù)[13]（autocorrelation function, ACF）和偏自相關(guān)系數(shù)[14]（partial autocorrelation function, PACF）.在參數(shù)估計(jì)階段，依據(jù)赤池信息準(zhǔn)則[15]（Akaike information criterion,AIC）和貝葉斯信息準(zhǔn)則[16]（Bayesian information criterion, BIC）確定最優(yōu)階數(shù).

1.2 多重共線(xiàn)性

多重共線(xiàn)性的產(chǎn)生是因?yàn)槟Ｐ椭凶宰兞恐g存在較高程度的相關(guān)性，所以自變量之間通?？梢韵嗷ゾ€(xiàn)性表出，即存在一組不全為0 的系數(shù)c0,c1,c2,··· ,cn，使得自變量在幾何平面上幾乎可以重疊成一條直線(xiàn)[17]，其計(jì)算公式為

多重共線(xiàn)性越嚴(yán)重，以最小二乘法估計(jì)的參數(shù)方差就會(huì)不斷增加，最終導(dǎo)致回歸方程的穩(wěn)健性[18]下降、個(gè)別自變量的顯著性關(guān)系不明顯等結(jié)果.目前，診斷多重共線(xiàn)性問(wèn)題方面最常用的指標(biāo)為容忍度（tolerance）T 和方差膨脹因子（variance inflation factor）Vi.兩者之間互為倒數(shù)，記xi的方差膨脹因子為Vi，其計(jì)算公式為

一旦診斷出模型存在多重共線(xiàn)性，就要想辦法解決該問(wèn)題.近幾年，很多學(xué)者對(duì)該領(lǐng)域進(jìn)行了大量的研究，在一定程度上緩解了共線(xiàn)性問(wèn)題帶來(lái)的預(yù)警準(zhǔn)確度不高的問(wèn)題.為進(jìn)一步提升預(yù)警準(zhǔn)確程度，本文選擇主成分分析法解決以上問(wèn)題.

1.3 主成分分析

主成分分析是將多種相關(guān)性較高的數(shù)據(jù)轉(zhuǎn)換為少數(shù)幾個(gè)綜合指標(biāo)，從而達(dá)到降維的目的.轉(zhuǎn)換后的指標(biāo)為原始指標(biāo)的線(xiàn)性組合，且轉(zhuǎn)換后的指標(biāo)數(shù)遠(yuǎn)少于轉(zhuǎn)換前的指標(biāo)數(shù)，此時(shí)也能保留原始數(shù)據(jù)的絕大部分信息.該過(guò)程[19]包括以下內(nèi)容：

假設(shè)原始數(shù)據(jù)包含的樣本量為n，每個(gè)樣本有p 個(gè)屬性，則原始數(shù)據(jù)構(gòu)成了一個(gè)n×p 維的矩陣，具體公式為

數(shù)據(jù)眾多屬性之間的單位和數(shù)量級(jí)通常存在差異，為了進(jìn)一步分析研究，應(yīng)先根據(jù)式(5)進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化變換

特征值越大對(duì)應(yīng)的方差貢獻(xiàn)率越大，于是通過(guò)累計(jì)方差貢獻(xiàn)率來(lái)確定最佳的主成分?jǐn)?shù)目.通常選取累計(jì)方差貢獻(xiàn)率達(dá)到80%的前n 個(gè)線(xiàn)性組合作為主成分[20]，其計(jì)算公式為

1.4 主成分回歸

假設(shè)隨機(jī)變量y 與主成分F1,F2,··· ,FP的線(xiàn)性回歸模型為[21]

式中，β0為回歸函數(shù)，常數(shù)項(xiàng)βp（i=1, 2,···, p）為自變量回歸系數(shù)，u 為隨機(jī)干擾項(xiàng).將式(8)進(jìn)行轉(zhuǎn)化后可得其矩陣形式為

根據(jù)最小二乘法對(duì)方程進(jìn)行偏導(dǎo)數(shù)極值求解，得到未知參數(shù)的最小二乘估計(jì)值為

2 算法模型構(gòu)建

在較大城市，移動(dòng)網(wǎng)絡(luò)目前幾乎實(shí)現(xiàn)了全網(wǎng)覆蓋，且手機(jī)的人群使用量也達(dá)到了較高水平，手機(jī)信令數(shù)據(jù)已經(jīng)成為城市客流聚集風(fēng)險(xiǎn)預(yù)測(cè)的一種常用數(shù)據(jù).運(yùn)營(yíng)商網(wǎng)絡(luò)具有以下兩點(diǎn)自身特性：1）運(yùn)營(yíng)商網(wǎng)絡(luò)覆蓋范圍和實(shí)際區(qū)域無(wú)法完全匹配，且誤差范圍不穩(wěn)定；2）某些區(qū)域存在基站同時(shí)覆蓋、信號(hào)漂移和干擾等因素.因此，運(yùn)營(yíng)商的手機(jī)信令數(shù)據(jù)往往與該區(qū)域的真實(shí)客流人數(shù)有較大差異，不能真實(shí)反映該區(qū)域的客流聚集程度.

本文從上海市內(nèi)客流密集的公共場(chǎng)所中選擇某區(qū)域作為研究對(duì)象，分析該區(qū)域運(yùn)營(yíng)商提供的手機(jī)用戶(hù)數(shù)據(jù)（該數(shù)據(jù)為經(jīng)過(guò)脫敏處理后的統(tǒng)計(jì)數(shù)據(jù)而非個(gè)體數(shù)據(jù).），選擇算法所需滯留用戶(hù)人數(shù)X1、滯留用戶(hù)7 天平均人數(shù)X2、新進(jìn)入用戶(hù)人數(shù)X3、新進(jìn)入用戶(hù)7 天平均人數(shù)X4、離開(kāi)用戶(hù)人數(shù)X5、離開(kāi)用戶(hù)7 天平均人數(shù)X6字段，對(duì)所選字段數(shù)據(jù)進(jìn)行預(yù)處理后得到用戶(hù)數(shù)分布及變化情況如圖1 所示.

圖1 手機(jī)用戶(hù)數(shù)分布圖Figure 1 Distribution of mobilephone users

2.1 ARIMA 模型參數(shù)選擇

模型預(yù)測(cè)前需要對(duì)輸入的時(shí)間序列進(jìn)行平穩(wěn)性檢驗(yàn).觀(guān)察圖1 可知該時(shí)間序列不平穩(wěn)，于是對(duì)時(shí)間序列一階進(jìn)行差分，得到自相關(guān)系數(shù)和偏自相關(guān)系數(shù)如圖2 和3 所示.

圖2 自相關(guān)系數(shù)分布圖Figure 2 Distribution of autocorrelation

圖3 偏自相關(guān)系數(shù)分布圖Figure 3 Distribution of partial autocorrelation

從圖2 和3 中可以看出自相關(guān)圖和偏自相關(guān)圖二階拖尾，依據(jù)AIC 和BIC 最小的原則可知ARIMA（2,1,2）為最優(yōu)模型.

模型參數(shù)確定后就可以檢驗(yàn)其殘差分布，合理的模型殘差分布應(yīng)該符合均數(shù)為0 的正太分布.本文繪制如圖4 所示的殘差分布圖，可見(jiàn)本文模型的殘差分布接近均值為0 的正態(tài)分布，符合作為對(duì)比實(shí)驗(yàn)的條件.

圖4 殘差分布圖Figure 4 Residual distribution

2.2 相關(guān)性分析與共線(xiàn)性診斷

分析運(yùn)營(yíng)商提供的手機(jī)用戶(hù)數(shù)據(jù)，得到數(shù)據(jù)之間的相關(guān)系數(shù)、容忍度和方差膨脹因子如表1 所示.相關(guān)系數(shù)矩陣表反映出數(shù)據(jù)之間的相關(guān)性非常高，除了X1和X2為0.88 之外，其余兩兩之間均超過(guò)了0.90，為防止后續(xù)建模中出現(xiàn)較為嚴(yán)重的多重共線(xiàn)性，需要診斷變量間的共線(xiàn)性.從表1 的計(jì)算結(jié)果中可以看出：除了變量X1和X2外，其他變量的容忍度T 均接近于0，方差膨脹因子則較大，X1的膨脹因子最小，但也超過(guò)了10.因此，根據(jù)容忍度[22]T 和方差膨脹因子[23]Vi這2 個(gè)指標(biāo)數(shù)值可以得出變量之間具有較嚴(yán)重共線(xiàn)性的結(jié)論.

表1 相關(guān)系數(shù)矩陣及共線(xiàn)性診斷Table 1 Correlation coefficient matrix and colinearity diagnosis

2.3 主成分回歸

本文在式（9）中加入結(jié)構(gòu)風(fēng)險(xiǎn)最小化參數(shù)，并結(jié)合噪聲和隨機(jī)參考項(xiàng)提出了以下訓(xùn)練模型：

該訓(xùn)練模型保留了主成分的回歸優(yōu)勢(shì)，同時(shí)簡(jiǎn)化了模型的復(fù)雜度，提升了模型的泛化率，進(jìn)一步加速了模型的運(yùn)行效率.

通過(guò)相關(guān)矩陣和特征值得出主成分碎石圖，如圖5 所示.從圖5 中可以看出：主成分1 和2 之間的拐角比較大，過(guò)渡非常陡峭，而其余主成分之間的過(guò)渡則比較平緩.從數(shù)值來(lái)看，相關(guān)矩陣的最大特征值約為6.0，第2 大特征值則小于1，而最小特征值接近于0，故從主成分碎石圖可以初步判定本文只需取一個(gè)主成分.

通過(guò)相關(guān)矩陣進(jìn)一步進(jìn)行主成分分析，得到如表2 所示的結(jié)果，第1 主成分的貢獻(xiàn)率已經(jīng)達(dá)到94.98%，第2 主成分的貢獻(xiàn)率為4.02%，第3 主成分的貢獻(xiàn)率已經(jīng)小于1%.從方差累計(jì)貢獻(xiàn)率來(lái)看，第1 主成分的累計(jì)貢獻(xiàn)率已經(jīng)達(dá)到了94.98%，即已經(jīng)可以解釋各類(lèi)原始數(shù)據(jù)總貢獻(xiàn)率的94.98%，故本文可以只取第1 主成分，因?yàn)榈? 主成分能反映原始變量接近95%的信息.主成分表達(dá)式為

本文先將原始數(shù)據(jù)進(jìn)行多重共線(xiàn)性診斷和主成分分析，再把提取出來(lái)的主成分與對(duì)應(yīng)時(shí)刻的真實(shí)客流數(shù)據(jù)進(jìn)行比對(duì)，通過(guò)逐步回歸擬合最終得出最優(yōu)的回歸表達(dá)式結(jié)果如表3 所示.調(diào)整后的可決系數(shù)高達(dá)99.15%，說(shuō)明整個(gè)模型擬合效果較好.主成分的一次項(xiàng)、二次項(xiàng)、三次項(xiàng)估計(jì)系數(shù)顯然都能通過(guò)方程參數(shù)t 檢驗(yàn)，說(shuō)明主成分的一次項(xiàng)、二次項(xiàng)、三次項(xiàng)均對(duì)因變量起重要作用.模型擬合的p 值遠(yuǎn)小于0.01，顯然也能通過(guò)F 檢驗(yàn)，故主成分回歸方程為

圖5 主成分碎石圖Figure 5 Principal component lithogram

表2 主成分分析Table 2 Principal component analysis

表3 主成分回歸結(jié)果Table 3 Result of principal component regression

根據(jù)表3 的檢驗(yàn)結(jié)果，預(yù)判定式（12）中最高階數(shù)為3 階，且3 階判定系數(shù)值高達(dá)99.2%，調(diào)整后的可決系數(shù)高達(dá)99.15%.判定系數(shù)和可決系數(shù)值較高，從側(cè)面反映了已定參數(shù)階數(shù)的可信度.又因?yàn)槟Ｐ蛿M合的p 值為2.2×10?6，遠(yuǎn)小于0.01，所以可確定主成分回歸方程定位3 階，系數(shù)值為表3 中Estimate 列對(duì)應(yīng)值，具體如式（13）所示.

3 預(yù)測(cè)效果評(píng)估及對(duì)比

3.1 實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)數(shù)據(jù)來(lái)源于運(yùn)營(yíng)商，數(shù)據(jù)采集開(kāi)始階段為2015年9 月至今，每條數(shù)據(jù)間隔10 min，數(shù)據(jù)總量截止到實(shí)驗(yàn)開(kāi)始階段為40×30×24×6 條.數(shù)據(jù)字段分別為區(qū)域名稱(chēng)（編號(hào)）、時(shí)間戳、用戶(hù)數(shù)（區(qū)域內(nèi)人數(shù)）、漫入數(shù)、漫出數(shù)5 個(gè)字段.本文截取了某區(qū)域的數(shù)據(jù)樣例展示，如表4 所示.

表4 中區(qū)域名是指每個(gè)區(qū)域的代號(hào)，時(shí)間指數(shù)據(jù)采集的具體時(shí)間點(diǎn)，用戶(hù)數(shù)為當(dāng)前基站監(jiān)測(cè)到的固定區(qū)域內(nèi)的人數(shù)，漫入數(shù)為當(dāng)前基站監(jiān)測(cè)到新加入用戶(hù)的數(shù)量，漫出數(shù)為當(dāng)前基站監(jiān)測(cè)到的從當(dāng)前基站流失人員的數(shù)量.

表4 某區(qū)域數(shù)據(jù)Table 4 Area data

3.2 參考指標(biāo)和誤差分布等級(jí)

將原始數(shù)據(jù)代入主成分回歸方程可得預(yù)測(cè)值，再將預(yù)測(cè)值與相應(yīng)真實(shí)值進(jìn)行對(duì)比，所得結(jié)果如圖6 和7 所示.預(yù)測(cè)值與真實(shí)值的最小誤差回歸PCA 為0.67%，平均誤差為25.1%，誤差中位數(shù)為7.79%，與ARIMA 模型相比平均誤差中位數(shù)明顯降低了.個(gè)別或少數(shù)奇異點(diǎn)的存在導(dǎo)致誤差平均值和誤差中位數(shù)值相差較大，其平均值不能充分反映誤差的整體情況，故應(yīng)以中位數(shù)為參考指標(biāo).

圖6 PCA 預(yù)測(cè)誤差分布Figure 6 Distribution of PCA prediction error

圖7 ARIMA 預(yù)測(cè)誤差分布Figure 7 Distribution of ARIMA prediction error

根據(jù)圖6 和7 將誤差范圍進(jìn)行分類(lèi)，定義誤差范圍低于10%的數(shù)量為1 級(jí)，定義誤差范圍在10%～20%之間的數(shù)量為2 級(jí)，定義誤差范圍在20%～30%之間的數(shù)量為3 級(jí)，定義誤差范圍高于30%的數(shù)量為4 級(jí).根據(jù)以上4 個(gè)等級(jí)劃分，針對(duì)回歸PCA 和ARIMA 分別畫(huà)出預(yù)測(cè)誤差分布圖，如圖8 所示.回歸PCA 的誤差范圍在10%以?xún)?nèi)的數(shù)量達(dá)到了進(jìn)行預(yù)測(cè)總數(shù)量的62%，誤差范圍超過(guò)30%的數(shù)量?jī)H為總數(shù)的9%；ARIMA 的誤差范圍在10%以?xún)?nèi)的數(shù)量?jī)H為34%，誤差超過(guò)30%的數(shù)量為15%，與回歸PCA 相比劣勢(shì)較為明顯，可見(jiàn)主成分回歸模型的整體預(yù)測(cè)效果較為理想.

3.3 實(shí)驗(yàn)對(duì)比

為說(shuō)明本文算法的可行性和高效性，本文設(shè)計(jì)了與ARIMA 算法的對(duì)比試驗(yàn)，對(duì)比效果如圖9 所示.

圖8 預(yù)測(cè)誤差級(jí)別分布Figure 8 Distribution of prediction error level

圖9 算法對(duì)比圖Figure 9 Algorithm comparison diagram

從圖9 中可以看出，本文算法在整體的預(yù)測(cè)準(zhǔn)確度方面與ARIMA 算法相比差距較小，在效率運(yùn)行方面如表5 所示.隨著迭代次數(shù)的增多，改進(jìn)PCA 算法的準(zhǔn)確度得到了提升，同時(shí)其運(yùn)行效率呈高斯分布，即呈現(xiàn)中間高兩頭低的趨勢(shì)，原因如下：隨著迭代次數(shù)的增多，模型復(fù)雜度逐步提升，模型訓(xùn)練效果更佳，泛化效果更好.然而，在復(fù)雜度提升的同時(shí)運(yùn)行效率略有降低，實(shí)驗(yàn)結(jié)果符合預(yù)期.

表5 運(yùn)行效率對(duì)比Table 5 Comparison of operating efficiency

4 結(jié) 語(yǔ)

本文從運(yùn)營(yíng)商提供的手機(jī)用戶(hù)數(shù)據(jù)出發(fā)，對(duì)每個(gè)維度潛在的特征和規(guī)律進(jìn)行分析和挖掘，發(fā)現(xiàn)各個(gè)維度之間存在較高共線(xiàn)性，若直接進(jìn)行回歸分析則會(huì)使預(yù)測(cè)結(jié)果受到質(zhì)疑.因此，借助主成分分析具有將多變量簡(jiǎn)化為少數(shù)不相關(guān)變量的優(yōu)點(diǎn)，有效克服了各個(gè)維度的共線(xiàn)性問(wèn)題，降低了算法的復(fù)雜度，同時(shí)與ARIMA 算法相比預(yù)測(cè)準(zhǔn)確度提升較為明顯.從該模型的分析結(jié)果中可以看出：因?yàn)槌浞挚紤]了運(yùn)營(yíng)商手機(jī)用戶(hù)數(shù)據(jù)的各個(gè)維度，所以模型的預(yù)測(cè)精度比較理想，可以為城市的有效管理、客流聚集風(fēng)險(xiǎn)的管控等提供有效的數(shù)據(jù)支撐.