李 威 魯鐵定 賀小星 錢文龍

1 東華理工大學測繪工程學院，南昌市廣蘭大道418號，330013 2 華東交通大學土木建筑學院，南昌市雙港東大街808號，330013 3 中鐵第一勘察設(shè)計院集團有限公司軌道交通工程信息化國家重點實驗室,西安市西影路2號,710043

近年來，不斷累積的IGS基準站坐標數(shù)據(jù)為大地測量學和地球動力學研究提供了豐富的數(shù)據(jù)基礎(chǔ)，也為GNSS坐標時間序列的研究提供了堅實的數(shù)據(jù)支撐[1]。GNSS坐標時間序列的精確預測對建筑物的變形監(jiān)測、氣候氣象的預測及區(qū)域內(nèi)地殼形變等研究具有重要意義。有研究指出，GNSS坐標時間序列在3個方向都有明顯的周期性運動，其中高程方向的季節(jié)性變化最為顯著[2-3]。另外，由于高程方向的噪聲分量往往大于水平方向，且噪聲模型的組成并不單一，使得建立高精度GNSS高程坐標時間序列預測模型的難度較大[4]。

諸多學者在GNSS坐標時間序列的擬合與預測領(lǐng)域進行了廣泛研究[5-8]，促進了GNSS坐標時間序列的研究進展，但針對高程坐標長期預測模型的研究仍存在一些不足：1)現(xiàn)階段時間序列模型多存在建模復雜、穩(wěn)定性差、無法動態(tài)調(diào)整模型參數(shù)等缺陷；2)模型對坐標數(shù)據(jù)的要求較高，建模精度容易受粗差影響，且未顧及數(shù)據(jù)缺失等不利影響；3)構(gòu)建的模型缺乏長期預測能力，在預測時需要滾動預測才能保持模型的精度，增加了模型的操作難度與運行時間。

針對這些問題，本文提出將Prophet模型應(yīng)用于GNSS高程坐標時間序列中，該模型對數(shù)據(jù)有良好的適應(yīng)能力，在長期預測中有較高的可靠性與精度，且無需太多關(guān)于時間序列的先驗知識。本文首先檢驗了Prophet模型在高程坐標時間序列中的適用性；再針對Prophet模型對非線性部分預測能力較弱的缺陷，引入隨機森林(RF)模型，得到Prophet-RF組合預測模型；最后根據(jù)2種模型的特性設(shè)計多種組合方案,并利用BJFS站的實測高程坐標數(shù)據(jù)檢驗不同組合方案的精度與適用性。

1 基本原理與評價指標

1.1 Prophet模型

Prophet模型由Taylor等[9]于2017年提出并發(fā)布了相關(guān)的開源軟件包，是一種可以分析時間序列的新模型，包括處理時間序列數(shù)據(jù)中的異常值和缺失值及對時間序列進行長短期預測并研究其變化規(guī)律。該模型在電離層異常探測[10]、市場銷量研究[11]、空氣質(zhì)量分析[12-13]等領(lǐng)域都有成功案例，但在GNSS坐標時間序列方面幾乎沒有涉及。

Prophet模型在預測時的循環(huán)分析共分為4個部分：建立模型、預測評估、表現(xiàn)問題、可視化檢測預測結(jié)果。使用Prophet模型分析數(shù)據(jù)時將4個部分構(gòu)成一個循環(huán)體，根據(jù)不同專業(yè)及數(shù)據(jù)的具體要求對模型進行人工調(diào)整，并結(jié)合模型的自動預測功能，使Prophet模型較傳統(tǒng)的時間序列模型有更好的靈活性與適用性[9]。

Prophet模型是一個廣義加法模型，并且使用基于貝葉斯的曲線擬合方法對時間序列進行平滑和預測[9-10]。模型由趨勢項、周期項及假日項組成，其加法模型的基本分解形式為：

y(t)=g(t)+s(t)+h(t)+ε(t)

(1)

其中，g(t)為趨勢項，表示時間序列中非周期部分的變化函數(shù)，包含不同程度的假設(shè)與可調(diào)節(jié)光滑度的參數(shù)[11]；s(t)為周期項，在Prophet模型中使用傅里葉級數(shù)對其進行構(gòu)造，常以周、月、季度或年作為單位；h(t)為假日項，表示因假期或特殊原因造成的不規(guī)則影響，在GNSS坐標時間序列中通常不存在此類情況，故在分解數(shù)據(jù)時可不必對假日項進行建模；ε(t)為殘差項，服從正態(tài)分布，表示模型中未預測到的隨機趨勢。具體描述見文獻[9]。

1.2 隨機森林模型

隨機森林(random forest，RF)模型是一種具備集成思想的機器學習模型，是由一定數(shù)量的弱學習器決策回歸樹(classification and regression tree，CART)結(jié)合bagging集成學習理論和隨機特征子空間構(gòu)成的強學習器[14-15]。RF模型的隨機性表現(xiàn)在隨機選擇樣本及隨機選擇特征，在選擇樣本時采用的是bootstrap重抽樣方法，其特點是既能保證每個樣本在每次抽取時都有相同的概率被抽中，又能在不降低樣本訓練規(guī)模的同時留出驗證集[16]；而在特征選擇上，CART使用基尼系數(shù)進行選擇，基尼系數(shù)越小表示特征越好。構(gòu)建RF模型的步驟如圖1所示。

圖1 隨機森林算法流程Fig.1 Random forest algorithm flow

RF模型的構(gòu)建流程為：

1)采用bootstrap重抽樣方法從原始數(shù)據(jù)集N中抽取ntree個訓練樣本集，每個訓練樣本集的大小約為原始數(shù)據(jù)集的三分之二。

2)對每個訓練樣本集分別建立CART，并從全部M個特征中隨機抽取mtry個特征作為CART的特征變量，再根據(jù)基尼系數(shù)最小原則選擇最佳方案對節(jié)點進行分裂。

3)每棵決策樹在生長期間都按照步驟2)進行，并保持mtry的大小不變；且每棵決策樹都要進行完全的生長，不對其進行任何剪枝操作。

4)由前3步生成的ntree棵決策樹組成一個隨機森林。在分類問題中，RF模型對檢驗集進行投票并根據(jù)少數(shù)服從多數(shù)的原則進行分類預測；而在回歸問題中，將每棵樹的結(jié)果等權(quán)相加，求得的均值作為最終預測值。

1.3 模型評價指標

本文使用均方根誤差(RMSE)與平均絕對誤差(MAE)作為各個模型擬合與預測精度的評判標準。RMSE與MAE的定義為：

(2)

(3)

式中，g為預測值，h為原始值。RMSE和MAE的值越小，表示模型的精度越高；反之，則表示模型的精度越低。

通過引入皮爾森相關(guān)系數(shù)(Pearson correlation coefficient)和斯皮爾曼相關(guān)系數(shù)(Spearman correlation coefficient)來檢驗各模型預測數(shù)據(jù)與原始數(shù)據(jù)的相關(guān)性。通常皮爾森相關(guān)系數(shù)對數(shù)據(jù)的要求較高，需要2個變量為線性關(guān)系且符合正態(tài)分布，而斯皮爾曼相關(guān)系數(shù)對數(shù)據(jù)沒有過多限制。根據(jù)表1可對數(shù)據(jù)的相關(guān)性強弱進行判斷[17]，相關(guān)系數(shù)的絕對值越接近1，表示其相關(guān)程度越強。

表1 相關(guān)性強弱判定表

2 Prophet-RF組合模型的構(gòu)建

選用BJFS站2008～2013年的高程坐標數(shù)據(jù)作為實驗數(shù)據(jù)，數(shù)據(jù)來源于中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(http://www.cgps.ac.cn/),采樣間隔為1/365.25 a，最大采樣頻率為365.25 Hz。本文所有實驗方案都保持相同的樣本集劃分：2008-01-01～2012-12-31的每日數(shù)據(jù)作為訓練樣本集，2013-01-01～12-31的數(shù)據(jù)作為測試樣本集。圖2為由Prophet模型得到的預測結(jié)果，可以看出，Prophet模型能夠較好地擬合原始數(shù)據(jù)中的周期性變化，且具有較好的抗粗差能力。

通過分析Prophet模型與RF模型的特性構(gòu)建不同的組合方案：方案1是基于Prophet模型自帶的分解功能，將數(shù)據(jù)分解為周期部分與非線性部分；方案2是假設(shè)在數(shù)據(jù)噪聲量級較大的情況下，使用Prophet模型對原始數(shù)據(jù)進行擬合；方案3是假設(shè)在數(shù)據(jù)噪聲量級較小的情況下，使用RF模型對原始數(shù)據(jù)進行擬合。

圖2 Prophet模型對BJFS站高程坐標的預測Fig.2 Prediction of elevation coordinate of BJFS station by Prophet model

2.1 Prophet-RF組合模型方案1

Prophet模型對線性趨勢有較好的擬合與預測效果，對于非線性部分的預測效果較差。根據(jù)這一特性，使用Prophet模型自帶的分解功能將高程坐標時間序列分解為趨勢項與周期項，其他部分看作殘差項。圖3中Prophet模型將BJFS站高程坐標時間序列分解為3個主要成分：1)由邏輯回歸函數(shù)或分段線性函數(shù)擬合的趨勢項；2)基于傅立葉級數(shù)構(gòu)建的以a為周期的季節(jié)項；3)由虛擬變量表示的以周為周期的季節(jié)項[13]。假日項可單獨添加，當數(shù)據(jù)中不存在此類情況時，采用默認分解即可。圖4為Prophet-RF組合模型方案1的流程。

圖3 Prophet模型分解坐標數(shù)據(jù)Fig.3 Prophet model decomposition coordinate data

圖4 Prophet-RF模型方案1流程Fig.4 Flow of scheme 1 of Prophet-RF model

Prophet-RF模型方案1的實現(xiàn)流程如下：

1)使用Prophet模型自帶的分解功能，將時間跨度為n天的坐標時間序列Sn分解為趨勢項trendn、以周為周期的季節(jié)項weeklyn、以a為周期的季節(jié)項yearlyn及殘差項errorn。

2)趨勢項、以周為周期的季節(jié)項及以a為周期的季節(jié)項分別使用Prophet模型建立預測模型，殘差項則使用RF模型建立預測模型。將每個分量模型的前t天作為訓練集，后k天作為測試集，且t+k=n。

3)將各分量模型得到的擬合結(jié)果與預測結(jié)果等權(quán)相加，得到Prophet-RF組合模型方案1的擬合值與預測值：

S擬合=trendt+weeklyt+yearlyt+errort

S預測=trendk+weeklyk+yearlyk+errork

(4)

2.2 Prophet-RF組合模型方案2

RF模型在某些噪聲較大的分類或回歸問題中可能會出現(xiàn)過擬合的情況，且數(shù)據(jù)中的噪聲量級往往難以判斷，故方案2在不明確數(shù)據(jù)噪聲量級的情況下，假設(shè)原始數(shù)據(jù)中噪聲量級較大。首先選用Prophet模型擬合坐標原始數(shù)據(jù)，再對擬合數(shù)據(jù)建立RF預測模型，避免了在數(shù)據(jù)噪聲量級較大時RF模型造成數(shù)據(jù)過擬合的缺點。圖5為Prophet-RF組合模型方案2的流程。

圖5 Prophet-RF模型方案2流程Fig.5 Flow of scheme 2 of Prophet-RF model

Prophet-RF模型方案2的實現(xiàn)流程如下：

1)將時間跨度為n天的原始坐標數(shù)據(jù)劃分為訓練樣本集和測試樣本集，前t天作為訓練樣本集，后k天作為測試樣本集，且t+k=n。

2)使用Prophet模型擬合訓練樣本集，并將得到的擬合數(shù)據(jù)看作新的特征變量。

3)利用RF模型對Prophet模型產(chǎn)生的特征變量建立預測模型，得到最終的擬合數(shù)據(jù)與預測數(shù)據(jù)。

2.3 Prophet-RF組合模型方案3

RF模型在數(shù)據(jù)噪聲含量較少的情況下能夠保持較好的獨立性，即不易過度擬合數(shù)據(jù)，故方案3在假設(shè)原始數(shù)據(jù)噪聲量級較小的情況下，使用RF模型擬合原始數(shù)據(jù)，避免了Prophet模型對非線性部分預測能力較弱的缺陷，從而能夠充分發(fā)揮Prophet模型的預測優(yōu)勢。圖6為Prophet-RF組合模型方案3的流程。

圖6 Prophet-RF模型方案3流程Fig.6 Flow of scheme 3 of Prophet-RF model

Prophet-RF模型方案3實現(xiàn)流程如下：

1)將時間跨度為n天的坐標原始數(shù)據(jù)劃分為訓練樣本集和測試樣本集，前t天作為訓練樣本集，后k天作為測試樣本集，且t+k=n。

2)使用RF模型擬合訓練樣本集，并將得到的擬合數(shù)據(jù)看作新的特征變量。

3)利用Prophet模型對RF模型產(chǎn)生的特征變量建立預測模型，并得到最終的擬合數(shù)據(jù)與預測數(shù)據(jù)。

2.4 方案對比與分析

2.4.1 精度對比

本文選用BJFS站高程坐標時間序列作為實驗數(shù)據(jù)，將Prophet模型及Prophet-RF組合模型方案1、方案2、方案3得到的擬合值與預測值分別與原始數(shù)據(jù)進行對比，同時引入經(jīng)典ARMA模型作為參照，得到各模型相應(yīng)的RMSE值與MAE值，結(jié)果見表2。

表2 不同模型的性能對比

由表2可知，ARMA模型的擬合精度最高，但預測精度最低，這是由于經(jīng)典ARMA模型容易出現(xiàn)過擬合現(xiàn)象, 且在長期預測中缺乏穩(wěn)定性。而單一的Prophet模型也能得到較高的擬合精度，但預測精度卻不如組合模型，表明模型的預測能力與擬合效果沒有直接關(guān)系。Prophet-RF組合模型方案1、方案2、方案3預測數(shù)據(jù)的RMSE值較單一的Prophet模型分別提升0.7%、14.7%、16.1%，MAE值分別提升7.8%、23.4%、24.2%，其中方案1的預測精度提升較小，而方案2及方案3的預測精度都有較大提升。

為檢驗各模型得到的擬合值、預測值與原始數(shù)據(jù)的關(guān)系，引入皮爾森相關(guān)系數(shù)和斯皮爾曼相關(guān)系數(shù)，由于ARMA模型在該數(shù)據(jù)中的適用性不強，故沒有分析其相關(guān)性。由表3可知，各模型的擬合值與原始數(shù)據(jù)具有強相關(guān)性，而預測值與原始數(shù)據(jù)具有中等強度相關(guān)性。通過對比分析表2和3的指標認為，組合模型方案3在BJFS站高程方向的預測精度最高、適用性最強。圖7為Prophet模型與Prophet-RF組合模型方案3預測誤差的對比，由預測值減去原始數(shù)據(jù)得到。

表3 不同模型的相關(guān)性對比

圖7 模型預測誤差對比Fig.7 Model prediction error comparison

2.4.2 組合模型分析

在方案1中，通過分解原始數(shù)據(jù)的組成成分，并根據(jù)Prophet模型和RF模型的特點，針對不同組成成分選用合適的預測模型，但該方案較單一的Prophet模型性能幾乎沒有提升。這可能是由于實測數(shù)據(jù)中往往存在由于各種誤差，如板塊運動、電離層延遲、觀測墩熱脹冷縮等造成的非線性變化，還有季節(jié)性海洋、陸地及大氣變化造成的地表水文負載等，這些都會對GNSS觀測數(shù)據(jù)產(chǎn)生負面影響[1]。所以方案1在分解數(shù)據(jù)時受到誤差的影響，無法較好地分離出各組成成分，從而導致模型精度提升較小。

方案2與方案3可看作一組對比實驗，分別在數(shù)據(jù)噪聲量級較大和較小的假設(shè)下進行。這2組方案較單一的Prophet模型預測精度都有較大提升，但方案3的預測精度更高，表明原始數(shù)據(jù)中的噪聲量級相對較小。方案2由于使用Prophet模型擬合原始數(shù)據(jù)時對非線性部分的擬合效果較差，再使用RF模型進行預測時精度可能會受到Prophet模型擬合效果的影響。而方案3使用RF模型擬合原始數(shù)據(jù)，避免了過擬合并對非線性部分進行了修正，再使用Prophet模型進行預測時能夠充分發(fā)揮模型優(yōu)勢，從而提高預測精度。

3 結(jié) 語

高精度的預測模型有助于研究GNSS坐標時間序列的長期變化與規(guī)律。本文將Prophet模型應(yīng)用于GNSS坐標數(shù)據(jù)中，該模型能夠分解坐標時間序列中的組成成分，且具有較強的可解釋性；在Prophet模型的基礎(chǔ)上引入RF模型，以解決Prophet模型對隨機部分預測能力較弱的缺點。實驗結(jié)果表明，Prophet-RF組合模型的預測精度相較于單一的Prophet模型，RMSE與MAE分別提升了16.1%與24.2%，并且對GNSS坐標時間序列有更好的長期預測能力。

由于數(shù)據(jù)往往受到噪聲及其他未知成分的影響，本文在進行組合模型擬合時充分考慮了數(shù)據(jù)中可能存在的情況，并結(jié)合各模型的優(yōu)缺點設(shè)計了多種組合方案，以保證實驗結(jié)果的可靠性，最終得到BJFS站高程方向的最優(yōu)組合方案。

本文對單個高程坐標時間序列進行了模型算法上的結(jié)合與研究，在今后的工作中應(yīng)當結(jié)合更多的模型或引入更新的算法，以進一步提高預測精度。