李 威 魯鐵定 賀小星 劉 瑞

1 東華理工大學(xué)測(cè)繪工程學(xué)院，南昌市廣蘭大道418號(hào)，330013 2 華東交通大學(xué)土木建筑學(xué)院，南昌市雙港東大街808號(hào)，330013

在GNSS坐標(biāo)數(shù)據(jù)處理中，數(shù)據(jù)缺失較為常見，不同的缺失情況對(duì)GNSS坐標(biāo)時(shí)間序列的分析會(huì)產(chǎn)生不同的影響[1]。拉格朗日方法適用于連續(xù)缺失數(shù)據(jù)在3個(gè)以內(nèi)的插值[2]，但隨著測(cè)量精度的提升，可能會(huì)出現(xiàn)插值效果不如簡(jiǎn)單線性插值的情況[3]。武艷強(qiáng)等[4]提出的三次樣條法能夠在一定程度上解決數(shù)據(jù)缺失較多的問題，但在數(shù)據(jù)連續(xù)缺失過多時(shí)，插值精度會(huì)大幅下降；邱榮海等[5]使用奇異譜迭代的插值方法提高了插值效率，但處理較為復(fù)雜，需要對(duì)嵌入窗口的插值階數(shù)進(jìn)行大量驗(yàn)證；尹玲等[6]使用神經(jīng)網(wǎng)絡(luò)進(jìn)行插值，該方法在觀測(cè)值缺失較多的情況下依然能夠得到較好的插值數(shù)據(jù)，但操作較為繁瑣；蘇利娜等[7]提出的基于模型和噪聲間相關(guān)性的插值方法對(duì)長(zhǎng)空缺數(shù)據(jù)有良好的插值效果，并發(fā)現(xiàn)數(shù)據(jù)空缺對(duì)周期性的影響較大，但該方法需要對(duì)時(shí)間序列具備大量的先驗(yàn)知識(shí)；蔡曉軍等[8]提出的多通道奇異譜插值方法不需要坐標(biāo)時(shí)間序列先驗(yàn)信息也能得到較高精度的插值數(shù)據(jù)，但交叉驗(yàn)證模型的參數(shù)且缺失數(shù)據(jù)需要用0填充可能帶來(lái)一定偏差。

本文針對(duì)GNSS坐標(biāo)時(shí)間序列中的插值問題，提出使用Prophet模型對(duì) GNSS坐標(biāo)時(shí)間序列進(jìn)行插值。該方法不需要關(guān)于時(shí)間序列的大量先驗(yàn)知識(shí)，并能在各種缺失情況下快速得到高精度的插值數(shù)據(jù)。

1 模型原理與數(shù)據(jù)處理

1.1 Prophet模型

Prophet模型插值的本質(zhì)是對(duì)缺失數(shù)據(jù)進(jìn)行預(yù)測(cè)，即在擬合的過程中自動(dòng)填補(bǔ)缺失值從而達(dá)到插值的效果[9]。Prophet模型克服了傳統(tǒng)插值方法的缺陷與限制，其在處理缺失值與異常值時(shí)魯棒性極強(qiáng)[10]。當(dāng)坐標(biāo)時(shí)間序列中存在缺失值時(shí)，Prophet模型能夠較好地確定GNSS坐標(biāo)時(shí)間序列的組成成分，并同時(shí)對(duì)具有多個(gè)季節(jié)性的周期數(shù)據(jù)進(jìn)行模擬[11]，使得到的插值數(shù)據(jù)能夠較好地契合時(shí)間序列的變化規(guī)律。

Prophet模型基于貝葉斯方法的曲線擬合來(lái)實(shí)現(xiàn)預(yù)測(cè)[11]，其可用廣義加法模型表示：

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

(1)

式中，g(t)為趨勢(shì)項(xiàng)，表示時(shí)間序列中非周期部分的變化趨勢(shì)；s(t)為周期項(xiàng)或季節(jié)項(xiàng)，通常以周或年為單位；h(t)為假日項(xiàng)，表示假日對(duì)1 d或多天內(nèi)不規(guī)則變化的影響；ε(t)為殘差項(xiàng)，表示模型未預(yù)測(cè)到的趨勢(shì)。

g(t)由邏輯回歸函數(shù)(logistic function)和分段線性函數(shù)(piecewise linear function)組成。由邏輯回歸函數(shù)表示g(t)為：

g(t)=C/(1+e-k(t-m))

(2)

式中，C為曲線的最大漸近值，即g(t)隨著t的增加趨于C；k為曲線的增長(zhǎng)率；m為曲線的偏移量。

s(t)使用傅里葉級(jí)數(shù)對(duì)其進(jìn)行構(gòu)造，表達(dá)形式為：

(3)

式中，T為期望時(shí)間序列具有的規(guī)則周期，當(dāng)T=365.25、N=10時(shí)，表示以年為周期；當(dāng)T=7、N=3時(shí)，表示以周為周期。

h(t)可用來(lái)反映時(shí)間序列中某時(shí)刻的特殊變動(dòng)，Prophet模型根據(jù)每個(gè)假日項(xiàng)在不同時(shí)刻下產(chǎn)生的影響構(gòu)建獨(dú)立的模型，并為各個(gè)假日項(xiàng)設(shè)置不同的前后窗口期，以及產(chǎn)生相應(yīng)的虛擬變量[10]。h(t)的表達(dá)形式為：

Z(t)=[1(t∈D1),…,1 (t∈DL)]

h(t)=Z(t)k,k～normal(0,γ)

(4)

式中，i為節(jié)假日，Ki為節(jié)假日的影響范圍；Di為i對(duì)應(yīng)的虛擬變量，表示第i個(gè)節(jié)假日的前后一段時(shí)間；L為時(shí)間序列中含有的節(jié)假日個(gè)數(shù)。

IGS基準(zhǔn)站的坐標(biāo)時(shí)間序列包含著測(cè)站的長(zhǎng)期線性變化趨勢(shì)，以及測(cè)站受地球物理效應(yīng)等外界因素影響造成的非線性變化，使測(cè)站位置產(chǎn)生周期性的震蕩變化[12-13]。所以當(dāng)坐標(biāo)時(shí)間序列中存在缺失值時(shí)，傳統(tǒng)的插值方法無(wú)法兼顧時(shí)間序列的周期變化。而Prophet模型能夠?qū)?shù)據(jù)中大量潛在的突變點(diǎn)進(jìn)行識(shí)別與監(jiān)測(cè)，再對(duì)趨勢(shì)變化的幅度進(jìn)行稀疏先驗(yàn)(與 L1 正則化效果相同)，并且Prophet模型采用的傅里葉級(jí)數(shù)可構(gòu)造適應(yīng)周期性變化的模型[14]，這使得Prophet模型在擬合時(shí)間序列的同時(shí)，能夠得到適應(yīng)時(shí)間序列周期變化的插值數(shù)據(jù)。圖1為Prophet模型對(duì)BJFS站坐標(biāo)數(shù)據(jù)的擬合圖，散點(diǎn)表示原始數(shù)據(jù)，曲線為擬合數(shù)據(jù)，陰影區(qū)域表示置信區(qū)間。由圖1可以看出，Prophet模型能夠較好地表現(xiàn)數(shù)據(jù)的變化趨勢(shì)以及坐標(biāo)時(shí)間序列的周期性，且沒有過度擬合原始數(shù)據(jù)，使得模型對(duì)粗差與離散點(diǎn)有較好的抵抗能力，這為得到高精度的插值數(shù)據(jù)提供了可能。

圖1 Prophet模型對(duì)BJFS站的數(shù)據(jù)擬合Fig.1 Data fitting of BJFS station using Prophet model

1.2 模擬數(shù)據(jù)的生成

幾乎所有的GNSS坐標(biāo)時(shí)間序列都表現(xiàn)出一個(gè)季節(jié)性的位移周期，它可以被模擬成一個(gè)周期項(xiàng)為1 a或0.5 a的4項(xiàng)傅里葉級(jí)數(shù)序列。最常見的軌跡模型是使用恒定的速度趨勢(shì)，在這種情況下，模擬數(shù)據(jù)遵循式(5)：

(5)

式中，xR為參考位置；tR為任意的參考時(shí)間，通常設(shè)置為平均觀測(cè)時(shí)間；v為測(cè)站的速度向量，可假定為常數(shù)；H為單位步長(zhǎng)函數(shù)；向量bj表示在時(shí)間tj上產(chǎn)生跳躍的方向和幅度，nJ為跳躍的次數(shù)；sk與ck為角頻率ωk的諧波傅里葉系數(shù)(每個(gè)位置向量對(duì)應(yīng)一個(gè))；nF為不同頻率的數(shù)量；角頻率ωk=2πτk，τk為角頻率對(duì)應(yīng)的周期。為了模擬年位移周期，選擇基本周期τ1=1 a，高次諧波周期τk=1/ka，這樣可確保由nF個(gè)正弦和nF個(gè)余弦構(gòu)成的循環(huán)每年僅重復(fù)一次。根據(jù)式(5)，分別模擬東、北、高程3個(gè)方向上的坐標(biāo)時(shí)間序列。通常GNSS坐標(biāo)時(shí)間序列的噪聲組成不能由單一的噪聲模型描述[15]，故模擬數(shù)據(jù)中包含線性趨勢(shì)+季節(jié)性信號(hào)，噪聲選用白噪聲與閃爍噪聲(WN+FN)的組合，且垂直分量大于水平分量。

2 Prophet模型插值實(shí)驗(yàn)

由于實(shí)測(cè)的GNSS坐標(biāo)數(shù)據(jù)包含各種誤差以及復(fù)雜的噪聲模型組合，需要通過模擬數(shù)據(jù)來(lái)檢驗(yàn)Prophet模型的插值能力。針對(duì)短周期和長(zhǎng)周期的坐標(biāo)時(shí)間序列，以及隨機(jī)缺失和連續(xù)缺失等情形，設(shè)計(jì)3組模擬實(shí)驗(yàn)。模擬實(shí)驗(yàn)1、2分別是在短周期和長(zhǎng)周期的坐標(biāo)時(shí)間序列中對(duì)隨機(jī)剔除的數(shù)據(jù)進(jìn)行插值，模擬實(shí)驗(yàn)3則對(duì)不同長(zhǎng)度的連續(xù)缺失數(shù)據(jù)進(jìn)行插值。在實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)中，分別進(jìn)行數(shù)據(jù)隨機(jī)缺失和連續(xù)缺失的插值實(shí)驗(yàn)。各實(shí)驗(yàn)均使用均方根誤差(RMSE)與平均絕對(duì)誤差(MAE)作為插值精度的評(píng)價(jià)指標(biāo)。

2.1 模擬實(shí)驗(yàn)1

模擬實(shí)驗(yàn)1主要針對(duì)短周期坐標(biāo)時(shí)間序列中隨機(jī)缺失的數(shù)據(jù)插值。根據(jù)式(5)得到模擬測(cè)站S1東、北、高程3個(gè)方向連續(xù)4 a共1 461 d的模擬數(shù)據(jù)。圖2為S1站3個(gè)方向的坐標(biāo)時(shí)間序列圖，將各方向時(shí)間序列隨機(jī)剔除5%、10%、15%的數(shù)據(jù)，再分別用3種方法對(duì)缺失數(shù)據(jù)進(jìn)行插值，并將得到的插值數(shù)據(jù)分別與缺失的原始數(shù)據(jù)進(jìn)行對(duì)比。

由表1中的精度指標(biāo)可以看出，3種方法在東向與北向的坐標(biāo)時(shí)間序列插值中都能保持較高的精度；而在高程方向上，盡管各方法的精度都有所下降，但仍能保持良好的插值效果。3種方法在3個(gè)方向上的插值精度并沒有隨著缺失比例的增加而出現(xiàn)明顯的降低，表明3種方法在短周期的時(shí)間序列中都有較好的插值效果。其中Prophet模型的插值精度最高，三次樣條法次之，拉格朗日法最差。

表1 S1站3種方法插值精度對(duì)比

2.2 模擬實(shí)驗(yàn)2

模擬實(shí)驗(yàn)2主要針對(duì)中長(zhǎng)周期坐標(biāo)時(shí)間序列隨機(jī)缺失的數(shù)據(jù)插值。由式(5)生成模擬測(cè)站S2東、北、高程3個(gè)方向連續(xù)10 a共3 653 d的模擬數(shù)據(jù)。圖3為S2站在3個(gè)方向上的時(shí)間序列圖，將各方向時(shí)間序列隨機(jī)剔除5%、10%、15%的數(shù)據(jù)，再分別采用拉格朗日法、三次樣條法與Prophet模型對(duì)缺失數(shù)據(jù)進(jìn)行插值，并將3種方法得到的插值數(shù)據(jù)與缺失的原始數(shù)據(jù)進(jìn)行對(duì)比。

圖3 S2站3個(gè)方向的坐標(biāo)時(shí)間序列Fig.3 Time series of three-dimensional coordinate at S2 station

從表2中的精度指標(biāo)可以看出，針對(duì)較長(zhǎng)周期的坐標(biāo)時(shí)間序列，3種模型都能較好地對(duì)缺失數(shù)據(jù)進(jìn)行插值。但拉格朗日法的插值精度波動(dòng)較大，這可能是由于當(dāng)缺失點(diǎn)較多時(shí)，拉格朗日插值多項(xiàng)式的次數(shù)增高，導(dǎo)致插值數(shù)據(jù)不穩(wěn)定；而Prophet模型與三次樣條法依然能夠提供高精度的插值數(shù)據(jù)。與模擬數(shù)據(jù)實(shí)驗(yàn)1的結(jié)果相同，Prophet模型在3種方法中插值精度最高。

表2 S2站3種方法插值精度對(duì)比

2.3 模擬實(shí)驗(yàn)3

模擬實(shí)驗(yàn)3主要檢驗(yàn)Prophet模型在數(shù)據(jù)連續(xù)缺失時(shí)的插值效果。由式(5)得到模擬測(cè)站S3東、北、高程3個(gè)方向連續(xù)6 a共2 192 d的模擬數(shù)據(jù)。圖4為S3站在3個(gè)方向上的時(shí)間序列圖，由于拉格朗日法對(duì)連續(xù)缺失數(shù)據(jù)的插值效果較差，故實(shí)驗(yàn)中不使用該方法。以2003-01-01為起點(diǎn)，分別剔除連續(xù)7 d、30 d、60 d、0.5 a、1 a、2 a共6種時(shí)間長(zhǎng)度的數(shù)據(jù)，再采用三次樣條法與Prophet模型對(duì)連續(xù)缺失數(shù)據(jù)進(jìn)行插值，并將2種方法得到的插值數(shù)據(jù)與缺失的原始數(shù)據(jù)進(jìn)行對(duì)比。由于三次樣條法在60 d時(shí)的插值精度已經(jīng)較低，故在缺失時(shí)間為0.5 a、1 a、2 a時(shí)只檢驗(yàn)Prophet模型的插值精度。

圖4 S3站3個(gè)方向的坐標(biāo)時(shí)間序列Fig.4 Time series of three-dimensional coordinate at S3 station

由表3的精度指標(biāo)可知，2種模型在東向與北向的插值精度都要優(yōu)于高程方向。其中，三次樣條法在數(shù)據(jù)連續(xù)缺失7 d和30 d時(shí)，能夠提供可用的插值數(shù)據(jù)；當(dāng)連續(xù)缺失數(shù)據(jù)為60 d時(shí)，其插值精度會(huì)大幅度下降。而Prophet模型在3種不同的缺失長(zhǎng)度下都能夠保持極高的插值精度；在缺失時(shí)間為0.5 a、1 a、2 a時(shí)，Prophet模型的插值精度并沒有隨著缺失時(shí)間的增加而明顯降低，仍能得到高精度的插值數(shù)據(jù)。

表3 S3站2種方法對(duì)連續(xù)缺失數(shù)據(jù)插值的精度

2.4 GNSS實(shí)測(cè)高程坐標(biāo)時(shí)間序列插值實(shí)驗(yàn)

為進(jìn)一步檢驗(yàn)Prophet模型插值方法的可靠性，將處理后得到的BJFS站(2008～2014年)高程方向坐標(biāo)時(shí)間序列(圖1)進(jìn)行插值研究與分析。數(shù)據(jù)來(lái)源于中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)(http://www.cgps.ac.cn/)，數(shù)據(jù)采樣間隔為1/365.25 a，采樣頻率為1/365.25 Hz。該站的時(shí)間跨度為6 a共2 192 d。按照5%、10%、15%的比例隨機(jī)剔除數(shù)據(jù)，然后對(duì)比Prophet模型、拉格朗日法以及三次樣條法對(duì)實(shí)測(cè)數(shù)據(jù)的插值效果。

由表4可知，3種插值方法在數(shù)據(jù)隨機(jī)缺失比例為5%、10%、15%時(shí)，對(duì)高程方向?qū)崪y(cè)數(shù)據(jù)的插值都能保持較高的精度，且插值效果十分穩(wěn)定。其中，Prophet模型的插值精度最高。

表4 BJFS站高程方向3種方法插值精度對(duì)比

當(dāng)BJFS站高程方向的實(shí)測(cè)數(shù)據(jù)連續(xù)缺失時(shí)，與模擬實(shí)驗(yàn)3相同，只選用Prophet模型與三次樣條法對(duì)缺失數(shù)據(jù)進(jìn)行插值。以2010-01-01為起點(diǎn)，分別剔除7 d、30 d、60 d、0.5 a、1 a、2 a共6種時(shí)間長(zhǎng)度數(shù)據(jù)，同樣在缺失時(shí)間為0.5 a、1 a、2 a時(shí)只檢驗(yàn)Prophet模型的精度(表5)。

表5 BJFS站高程方向2種方法對(duì)連續(xù)缺失數(shù)據(jù)插值的精度

由表5可知，Prophet模型與三次樣條法在實(shí)測(cè)數(shù)據(jù)連續(xù)缺失7 d時(shí)，都能夠得到較高的插值精度；而當(dāng)實(shí)測(cè)數(shù)據(jù)連續(xù)缺失30 d和60 d時(shí)，2種方法的插值精度都有不同程度的下降，但Prophet模型的插值精度遠(yuǎn)高于三次樣條法；在連續(xù)缺失時(shí)間為0.5 a、1 a、2 a時(shí)，Prophet模型的精度只有小幅度下降，且表現(xiàn)得十分穩(wěn)定。

圖5是Prophet模型與三次樣條法對(duì)BJFS站高程方向上連續(xù)缺失數(shù)據(jù)的插值對(duì)比圖?？梢钥闯?，三次樣條法在連續(xù)空缺較多時(shí)插值數(shù)據(jù)會(huì)偏離原始數(shù)據(jù)，從而產(chǎn)生較大的偏差；而Prophet模型得到的插值數(shù)據(jù)更加符合原始數(shù)據(jù)的變化。圖6(2010年)是使用Prophet模型在連續(xù)缺失0.5 a、1 a、2 a時(shí)的插值圖?？梢钥闯?，在缺失量較大時(shí)，Prophet模型的插值效果依然穩(wěn)定，并能夠較好地體現(xiàn)原始數(shù)據(jù)的趨勢(shì)變化與周期性。

圖5 BJFS站不同方法的插值對(duì)比Fig.5 Interpolation comparison of different methods at BJFS station

圖6 BJFS站高程方向的Prophet模型插值Fig.6 Prophet model interpolation of elevation direction at BJFS station

2.5 Prophet模型的靈活度

使用Prophet模型插值時(shí)，時(shí)間序列趨勢(shì)的變化靈活性不足(即擬合不足)或過于靈活(即過度擬合)，導(dǎo)致模型的默認(rèn)參數(shù)(默認(rèn)靈活度為0.05)有時(shí)不能得到最優(yōu)的插值數(shù)據(jù)，這時(shí)需要手動(dòng)調(diào)整模型的靈活性以調(diào)整稀疏先驗(yàn)的程度。

以模擬測(cè)站S1與BJFS站高程方向坐標(biāo)時(shí)間序列為例分析靈活度對(duì)Prophet模型插值的影響。按照15%的比例隨機(jī)剔除數(shù)據(jù)，再選擇常用的靈活度參數(shù)得到Prophet模型在不同靈活度下的插值數(shù)據(jù)，最后將每次得到的插值數(shù)據(jù)與原始數(shù)據(jù)對(duì)比得到相應(yīng)的RMSE。圖7中，x軸表示選取的靈活度參數(shù)，y軸為對(duì)應(yīng)的RMSE?？梢钥闯?，通常在靈活度為0.001時(shí)模型的插值效果最差，當(dāng)靈活度上升到0.05時(shí)插值效果會(huì)有較大提升。插值精度往往會(huì)隨著靈活度的上升而提高；但達(dá)到一定的精度后，即使靈活度上升，精度也不再提高，而是逐漸趨于平穩(wěn)甚至有小幅下降。

圖7 Prophet模型靈活度對(duì)精度的影響Fig.7 The influence of the flexibility of the Prophet model on accuracy

3 結(jié) 語(yǔ)

1)在建模過程中，由于Prophet模型不會(huì)過度擬合原始數(shù)據(jù)，所以能夠降低粗差或離群值對(duì)模型的影響。Prophet模型在GNSS坐標(biāo)時(shí)間序列中具有較好的適用性，這為Prophet模型的高精度插值提供了可能性。

2)Prophet不是單純的數(shù)學(xué)模型，它能夠兼顧數(shù)據(jù)的趨勢(shì)變化，從而提供良好的插值效果。Prophet模型在短周期與中長(zhǎng)周期的坐標(biāo)時(shí)間序列中都能保持良好的插值效果。而當(dāng)數(shù)據(jù)連續(xù)缺失時(shí)，Prophet模型較傳統(tǒng)模型優(yōu)勢(shì)更明顯。在小比例和非連續(xù)缺失的情況下，仍可選用拉格朗日法與三次樣條法進(jìn)行插值，但Prophet模型通常能夠提供更高精度的插值數(shù)據(jù)。

3)Prophet模型的優(yōu)勢(shì)在于簡(jiǎn)便性與快捷性，其不需要太多關(guān)于時(shí)間序列的先驗(yàn)知識(shí)，盡管它可以根據(jù)專業(yè)需求來(lái)調(diào)整模型以達(dá)到更好的效果，但在默認(rèn)設(shè)置下仍能給予較好的反饋。所以處理數(shù)據(jù)時(shí)不用過多地關(guān)注模型的選擇，而可以將重點(diǎn)放在對(duì)數(shù)據(jù)的分析中。

4)當(dāng)對(duì)數(shù)據(jù)插值的精度要求較高時(shí)，Prophet模型有時(shí)不能一次得出理想的結(jié)果，這時(shí)需要不斷調(diào)整模型的靈活度，而模型本身又無(wú)法自動(dòng)給出最優(yōu)解，所以盡管可以根據(jù)經(jīng)驗(yàn)較快地鎖定最優(yōu)靈活度的范圍，但仍會(huì)增加數(shù)據(jù)處理的時(shí)間。