基于時(shí)間序列模型的青藏鐵路路基變形預(yù)測(cè)

陳衛(wèi)雄

(中鐵十二局集團(tuán)鐵路養(yǎng)護(hù)工程有限公司, 拉薩 250014)

青藏鐵路沿線經(jīng)過(guò)的多年凍土區(qū)地質(zhì)條件十分復(fù)雜,不良凍土現(xiàn)象發(fā)育,至青藏鐵路2001年修建時(shí)，穿越多年凍土區(qū)約為546.4 km,其中高溫極不穩(wěn)定區(qū)和高溫不穩(wěn)定區(qū)長(zhǎng)274.25 km，低溫基本穩(wěn)定區(qū)和低溫穩(wěn)定區(qū)長(zhǎng)170.48 km，融區(qū)長(zhǎng)101.68 km。在所有凍土區(qū)線路中，高含冰量地段223.16 km，低含冰量地段221.57 km。

在青藏高原多年凍土區(qū)，多年凍土表對(duì)巖土工程的穩(wěn)定性有顯著影響。熱邊界和土壤性質(zhì)是影響多年凍土表的關(guān)鍵因素。復(fù)雜的地質(zhì)環(huán)境和人類活動(dòng)會(huì)導(dǎo)致熱邊界和土壤性質(zhì)的不確定性[1]。凍土對(duì)路基高程的影響主要體現(xiàn)在，由于溫度和承載力的變化，達(dá)到了路基強(qiáng)度、變形的協(xié)調(diào)條件，本質(zhì)上是水、熱、力、變形的耦合作用。其對(duì)路基變形的形象，大大威脅了路基的安全性與穩(wěn)定性。青藏鐵路沿線已出現(xiàn)如K1401+888橋7號(hào)墩梁縫擠死、K1159+833格臺(tái)限位塊開(kāi)裂、K1260+821第5片梁格端抗震樁拔起、K0973+551格臺(tái)墊石裂縫等危險(xiǎn)狀況。

由于凍土問(wèn)題的復(fù)雜性，帶來(lái)問(wèn)題的嚴(yán)重性，近年有諸多學(xué)者展開(kāi)了對(duì)鐵路凍土融沉、凍脹等的預(yù)測(cè)研究[2-9]。其中，苗姜龍[3]與祁長(zhǎng)青等[4]基于BP(back propagation)神經(jīng)網(wǎng)絡(luò)方法，建立了多變形因素的預(yù)測(cè)模型，并采用COMSOL數(shù)值模擬軟件進(jìn)行二次開(kāi)發(fā)。楊笑男[5]采用層次分析的方法，依托GIS(geographic information system)技術(shù)進(jìn)行地質(zhì)災(zāi)害危險(xiǎn)性的評(píng)級(jí)與預(yù)測(cè)。Chen等[8]引入了人工神經(jīng)網(wǎng)絡(luò)(artificial neural network,ANN)和長(zhǎng)短期記憶((long short-term memory, LSTM)網(wǎng)絡(luò)兩種基于深度學(xué)習(xí)技術(shù)的模型來(lái)預(yù)測(cè)鐵路路基凍脹變形。由于隨時(shí)間變化，溫度隨之改變，將會(huì)對(duì)路基造成影響，王文等[10]基于事件概率回歸估計(jì)，對(duì)青藏鐵路沿線平均最高最低氣溫變化趨勢(shì)進(jìn)行了預(yù)測(cè)與概率估計(jì)。李棟梁等[11]預(yù)測(cè)了青藏高原及鐵路沿線未來(lái)50年平均最低、最高和年地表溫度。此外，常斌等[12]應(yīng)用神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型對(duì)青藏鐵路沿線凍土通風(fēng)管路基溫度場(chǎng)進(jìn)行了分析。Zhang等[13]提出了基于季凍區(qū)3年現(xiàn)場(chǎng)監(jiān)測(cè)溫度數(shù)據(jù)的經(jīng)驗(yàn)?zāi)Ｐ?，表明路基溫度分布主要受地理位置和路基高度影響，橫向溫差普遍存在。對(duì)于養(yǎng)護(hù)措施，Liu等[14]研究表明碎石夾層路堤與垂直通風(fēng)管道相結(jié)合可以有效控制多年凍土區(qū)地溫，提高高速公路路堤穩(wěn)定性。周亞龍等[15]基于三維有限元分析，預(yù)測(cè)了輸電塔熱棒樁基的長(zhǎng)期降溫效果，并給出了工程建議。Ni等[16]使用沉降指數(shù)、風(fēng)險(xiǎn)區(qū)帶指數(shù)和許用承載力指數(shù)來(lái)評(píng)估青藏高原多年凍土區(qū)的沉降風(fēng)險(xiǎn)。Duvillard等[17]評(píng)估法國(guó)阿爾卑斯山永久凍土上基礎(chǔ)設(shè)施的不穩(wěn)定風(fēng)險(xiǎn)水平，并討論利益相關(guān)者制定的適應(yīng)和緩解策略。

現(xiàn)直接從數(shù)據(jù)角度出發(fā)，結(jié)合統(tǒng)計(jì)推斷，通過(guò)時(shí)間序列模型對(duì)青藏鐵路路基高程數(shù)據(jù)進(jìn)行分析，由Box-Jenkins建模方法確定合適的時(shí)間序列模型，并由此進(jìn)行青藏鐵路路基高程至2023年危險(xiǎn)點(diǎn)的預(yù)測(cè)，并由預(yù)測(cè)結(jié)果給出工程建議。

1 青藏鐵路路基高程數(shù)據(jù)預(yù)處理

K1425+050處左側(cè)各年月數(shù)據(jù)為例,如表1所示。

表1 K1425+050處左側(cè)各月路基高程數(shù)據(jù)

由于青藏鐵路路基位處凍土區(qū)區(qū)段部分受惡劣環(huán)境影響，從匯總各年總?cè)诔亮康臄?shù)據(jù)來(lái)看，各年的路基高程呈下降趨勢(shì)，因此數(shù)據(jù)有明顯的趨勢(shì)性，下面通過(guò)差分法消除數(shù)據(jù)趨勢(shì)性，為下一步模型建立提供基礎(chǔ)。

選取差分法去除路基高程數(shù)據(jù)的趨勢(shì)性，設(shè)時(shí)間序列{Yt}的觀察序列為{yt,t=1,2,…,N}，則序列的一階差分即為

?Yt=Yt-Yt-1=(1-B)Yt

(1)

式(1)中：?稱為差分算子；B為后移算子，即BYt=Yt-1。差分后的數(shù)據(jù)需經(jīng)過(guò)DF檢驗(yàn)(Dickey-Fuller test)確定其是否平穩(wěn)，即檢驗(yàn)γ是否為零。設(shè)差分后DF檢驗(yàn)?zāi)Ｐ蜑?/p>

?Yt=a1+a2t+γYt-1+ut

(2)

式(2)中：a1、a2為外生變量；t為時(shí)間；ut為白噪聲序列。對(duì)差分后數(shù)據(jù)進(jìn)行統(tǒng)計(jì)特征分析，畫出一階差分直方圖，如圖1所示。

圖1 一階差分?jǐn)?shù)據(jù)直方圖Fig.1 Histogram of first-order difference data

經(jīng)計(jì)算可得，一階差分后，K1425+050處左側(cè)數(shù)據(jù)平均值為-0.001 233 645，標(biāo)準(zhǔn)差為0.003 780 734，樣本偏度為0.878 933 3，樣本峰度為7.432 754，樣本均值在(0±2) s，即2倍標(biāo)準(zhǔn)差范圍內(nèi)，可以認(rèn)為一階差分后的數(shù)據(jù)為零均值過(guò)程，符合Box-Jenkins方法。

2 青藏鐵路路基高程數(shù)據(jù)時(shí)間序列模型識(shí)別與初步定階

路基高程一階差分?jǐn)?shù)據(jù)的自相關(guān)系數(shù)圖與偏相關(guān)系數(shù)圖，如圖2與圖3所示。

圖2 一階差分自相關(guān)系數(shù)圖Fig.2 The first-order differential autocorrelation coefficient diagram

圖3 一階差分偏相關(guān)系數(shù)圖Fig.3 The first-order differential partial correlation coefficient diagram

(3)

顯然有

所以只需檢驗(yàn)

i=1,2,…,M-q

(4)

是否成立，即可判斷零均值平穩(wěn)序列{Xt}的自相關(guān)系數(shù)是否在q步截尾，是否建立MA模型(moving average model)，其中M是自相關(guān)階數(shù)。

因此，可認(rèn)為自相關(guān)系數(shù)2階截尾，可考慮MA(2)建模。

(5)

所以只需檢驗(yàn)

(6)

是否成立，即可判斷零均值平穩(wěn)序列{Xt}自相關(guān)系數(shù)是否在p步截尾，是否建立AR模型(auto regressive model)，M為偏相關(guān)系數(shù)。觀察圖3可知各偏相關(guān)系數(shù)滿足

(7)

因此，可選擇建立AR(1)模型。

通過(guò)Box-Jenkins方法建立進(jìn)一步的ARMA模型，其特點(diǎn)為初步設(shè)定相關(guān)模型為ARMA(n,n-1)，即初步設(shè)定的相關(guān)模型中自回歸的階數(shù)比移動(dòng)平均的階數(shù)高一階，而選用ARMA(n,n-1)模型來(lái)擬合時(shí)間序列的原因如下[18]。

(1)AR(p)、MA(q)、ARMA(p,q)模型均為ARMA(n,n-1)模型的特殊情形。

(2)對(duì)任一平穩(wěn)隨機(jī)系統(tǒng)，均可用ARMA(n,n-1)模型近似。

(3)以差分方程理論可證明，若自回歸的階數(shù)是n，則移動(dòng)平均的階數(shù)是n-1。

(4)從連續(xù)系統(tǒng)的離散化過(guò)程來(lái)看，ARMA(n,n-1)具有合理性。在一個(gè)自回歸階數(shù)為n階，移動(dòng)平均階數(shù)為任意階的線性微分方程的形式下，如果對(duì)一個(gè)連續(xù)自回歸移動(dòng)平均過(guò)程進(jìn)行一致區(qū)間上的抽樣，那么，抽樣過(guò)程的結(jié)果是ARMA(n,n-1)。

用Box-Jenkins方法建立模型，即從低階到高階對(duì)模型進(jìn)行擬合及檢驗(yàn)，若{Xt}屬于ARMA模型，則采用ARMA(n,n-1)，并從n=2到高階對(duì)模型逐個(gè)進(jìn)行擬合和檢驗(yàn)，根據(jù)AIC準(zhǔn)則最終確定模型。對(duì)選出的模型進(jìn)行適應(yīng)性檢驗(yàn)，進(jìn)一步檢驗(yàn)?zāi)Ｐ汀?/p>

結(jié)合上述模型思想，可以嘗試建立ARMA(2,1)模型進(jìn)行擬合。下面對(duì)選出的3個(gè)模型進(jìn)行參數(shù)的擬合和模型選擇。

3 基于AIC準(zhǔn)則的模型確定方法

對(duì)于平穩(wěn)序列{Xt,t=0,±1,±2,…}，通過(guò)R軟件用極大似然估計(jì)可算得選出的ARIMA(0,1,1)、ARIMA(1,1,0)、ARIMA(2,1,1)模型參數(shù)值，且各模型的AIC值也可計(jì)算得到，結(jié)果如表2所示。

表2 模型參數(shù)估計(jì)值與各AIC值

由表2可知，最小AIC值在模型為ARIMA(2,1,1)時(shí)取得，因此，確定最終模型為ARIMA(2,1,1)。模型表達(dá)為

(8)

ARIMA模型是分析和研究時(shí)間序列問(wèn)題的常用方法。該方法能更有效地表征數(shù)據(jù)序列的線性特性，具有較好的短期預(yù)測(cè)效果[19]。

4 模型的適應(yīng)性檢驗(yàn)與預(yù)測(cè)

從殘差來(lái)看，時(shí)間序列模型得到的殘差應(yīng)具有純隨機(jī)型，因此，需要判斷殘差的正態(tài)性，通過(guò)圖4所示殘差數(shù)據(jù)的Q-Q圖(quantile-quantile plot)可以看出，殘差散點(diǎn)集中于直線附近，因此模型(4)對(duì)應(yīng)的殘差具備正態(tài)性。

圖4 殘差數(shù)據(jù)的Q-Q圖Fig.4 Q-Q diagram of residual data

通過(guò)自相關(guān)系數(shù)圖與偏相關(guān)系數(shù)圖，如圖5、圖6所示，與式(2)與式(3)判斷，殘差序列為純隨機(jī)序列。因此，模型(4)具備適應(yīng)性。因此由該模型進(jìn)行預(yù)測(cè)，預(yù)測(cè)結(jié)果如圖7所示。

圖5 殘差自相關(guān)系數(shù)圖Fig.5 Residual autocorrelation coefficient plot

由圖7可見(jiàn)，模型對(duì)實(shí)際數(shù)據(jù)的擬合較好，而預(yù)測(cè)數(shù)據(jù)大致呈現(xiàn)直線狀態(tài)，這是由于所選測(cè)點(diǎn)處，數(shù)據(jù)線性趨勢(shì)較為明顯。

圖7 路基高程預(yù)測(cè)圖Fig.7 Roadbed elevation prediction map

預(yù)測(cè)數(shù)據(jù)與測(cè)試集數(shù)據(jù)記錄于表3中。

由表3可見(jiàn)，由于數(shù)據(jù)具有明顯的線性性質(zhì)，預(yù)測(cè)值呈現(xiàn)直線下降趨勢(shì)，無(wú)法真實(shí)地模擬數(shù)據(jù)的情況，這是由于青藏鐵路氣候環(huán)境復(fù)雜，影響其高程變化的因素較多，不會(huì)呈現(xiàn)嚴(yán)格的統(tǒng)計(jì)性質(zhì)，且從預(yù)測(cè)絕對(duì)值來(lái)看，預(yù)測(cè)值沒(méi)有呈現(xiàn)偏差越來(lái)越大的趨勢(shì)，因此，預(yù)測(cè)結(jié)果具備可信度。

表3 路基高程預(yù)測(cè)值與真實(shí)值比較

按照以上建模步驟建立全部測(cè)點(diǎn)的融沉值的時(shí)間序列模型，并進(jìn)行預(yù)測(cè)，具體步驟如下。

Step 1由于路基融沉值在不同測(cè)點(diǎn)處的數(shù)據(jù)總是呈現(xiàn)一定趨勢(shì)性的，因此默認(rèn)進(jìn)行一階差分來(lái)消除數(shù)據(jù)趨勢(shì)性。

Step 2當(dāng)差分后數(shù)據(jù)仍不是平穩(wěn)數(shù)據(jù)時(shí)，為判斷數(shù)據(jù)是否具有周期性，從3階差分逐步進(jìn)行到12階差分，每次差分后判斷序列是否平穩(wěn)，以此找出序列周期性。

Step 3找到序列周期性后，對(duì)數(shù)據(jù)進(jìn)行去除周期性以使數(shù)據(jù)平穩(wěn)，對(duì)平穩(wěn)化后的數(shù)據(jù)進(jìn)行模型定階，并以此模型進(jìn)行路基高程數(shù)據(jù)的預(yù)測(cè)。

結(jié)合上述步驟得到各路基測(cè)點(diǎn)處路基高程至2023年12月預(yù)測(cè)值，根據(jù)其沉降量與路基左右側(cè)差值，總結(jié)了各預(yù)測(cè)中可能出現(xiàn)重大變形的危險(xiǎn)點(diǎn)如表4所示，在預(yù)測(cè)期內(nèi)，測(cè)點(diǎn)左右兩側(cè)路基高程差值出現(xiàn)較大差值的測(cè)點(diǎn)，如表4所示。為工程維護(hù)等提供參考。

表4中為各測(cè)點(diǎn)路基變形值，路基左右側(cè)差值在預(yù)測(cè)期內(nèi)前十的測(cè)點(diǎn)，由兩側(cè)差值較大的測(cè)點(diǎn)中可見(jiàn)，集中于K1476+550、K1476+650、K1476+600路基的左右兩側(cè)出現(xiàn)明顯長(zhǎng)距離的差異，因此，在K1476+600地點(diǎn)附近路基需要得到重點(diǎn)維護(hù)。

表4 預(yù)測(cè)期出現(xiàn)危險(xiǎn)變形測(cè)點(diǎn)

5 結(jié)論

青藏高原獨(dú)特的地理環(huán)境導(dǎo)致多年凍土區(qū)沿線的路基都不同程度地出現(xiàn)了凍脹、融沉等地質(zhì)災(zāi)害，所以研究多年凍土區(qū)路基沉降變形十分有必要。以青藏鐵路多年凍土區(qū)路基變形的數(shù)據(jù)為研究對(duì)象，得到如下結(jié)果。

(1)Box-Jenkins建模方法確定合適的時(shí)間序列模型為ARIMA(2,1,1)，通過(guò)了模型適應(yīng)性檢驗(yàn)，模型對(duì)實(shí)際數(shù)據(jù)擬合較好，預(yù)測(cè)結(jié)果具備可信度。

(2)總結(jié)了預(yù)測(cè)中可能出現(xiàn)重大變形的10個(gè)危險(xiǎn)點(diǎn)，在預(yù)測(cè)期內(nèi)，測(cè)點(diǎn)左右兩側(cè)路基高程差值出現(xiàn)較大差值的10個(gè)測(cè)點(diǎn)，為工程維護(hù)等提供參考。

(3)測(cè)點(diǎn)K1476+600附近，路基的左右兩側(cè)出現(xiàn)明顯長(zhǎng)距離的差異，該地點(diǎn)需要得到重點(diǎn)維護(hù)。

研究結(jié)果可為工程設(shè)計(jì)、施工和維護(hù)提供有價(jià)值的參考，并為青藏高原多年凍土融化沉降危害的預(yù)警和預(yù)防提供參考。但針對(duì)長(zhǎng)期預(yù)測(cè)效果精度的研究，還需進(jìn)一步擴(kuò)展。