李世璽，孫憲坤，尹 玲，張仕森

0 引言

利用全球定位系統(tǒng)（global position system，GPS）連續(xù)觀測高程數(shù)據(jù)進(jìn)行噪聲分析或形變監(jiān)測時，缺失數(shù)據(jù)的數(shù)量以及缺失數(shù)據(jù)的連續(xù)程度都會對分析

結(jié)果產(chǎn)生較大影響[1]，因此需要提高缺失數(shù)據(jù)的補(bǔ)全或預(yù)測精度，滿足分析需求，從而保證分析的可靠性?；谏窠?jīng)網(wǎng)絡(luò)的預(yù)測方法與一些傳統(tǒng)方法相比，則具有很強(qiáng)的非線性特征學(xué)習(xí)能力[2-5]。如文獻(xiàn)[6]利用長短期記憶神經(jīng)網(wǎng)絡(luò)（long short-term memory，LSTM）長時間的特征記憶能力對上海佘山站的高程坐標(biāo)進(jìn)行了預(yù)測，其30 d滾動預(yù)測結(jié)果的均方根誤差約為 2.4 mm，擬合趨勢也較為準(zhǔn)確，但使用了缺失數(shù)據(jù) 2端所有的可用數(shù)據(jù)訓(xùn)練模型，不符合實(shí)際應(yīng)用中的預(yù)測要求，且缺乏對GPS時間序列自身特性的分析。而在時間序列預(yù)測問題中，采用結(jié)合了時間序列數(shù)據(jù)自身特點(diǎn)的混合預(yù)測模型，往往能夠取得更好的預(yù)測效果，而非某種預(yù)測模型的簡單應(yīng)用[7-10]。

由于GPS高程時間序列是1種非平穩(wěn)非線性的1維時間序列，它包含了確定的長期趨勢、周期性的震蕩模式，以及短期的非線性變化規(guī)律和噪聲[11]，其中在我國境內(nèi)噪聲通常為白噪聲（white noise，WN）加閃爍噪聲（flicker noise，F(xiàn)N）[12-13]。遍布于整個時間序列的噪聲，使模型的整體預(yù)測結(jié)果具有不確定性，不同頻率的信號相互交疊既掩蓋了序列原本的變化規(guī)律，同時也增加了重復(fù)的時間序列信息，它們都使模型學(xué)習(xí)序列特征更加困難。因此提高時間序列的預(yù)測精度，應(yīng)先將噪聲分離，并從原始的1維序列中挖掘更多有效的特征模式用于預(yù)測。

混沌時間序列預(yù)測方法是著眼于分析時間序列自身特點(diǎn)的 1類非線性時間序列預(yù)測方法，它首先通過降噪減少序列的無用信息，再利用序列各時間節(jié)點(diǎn)間的相關(guān)性以及序列片段包含信息量的多少，將1維時間序列映射到高維相空間，從而恢復(fù)出時間序列系統(tǒng)特征，最后使用多特征序列進(jìn)行預(yù)測。混沌預(yù)測方法在許多領(lǐng)域都已成功應(yīng)用，如文獻(xiàn)[14]將 1維的徑流量時間序列進(jìn)行相空間重構(gòu)，恢復(fù)了復(fù)雜水文系統(tǒng)的時空特征，并將重構(gòu)后的多維序列作為預(yù)測模型的輸入進(jìn)行訓(xùn)練，提高了長江徑流量預(yù)測的穩(wěn)定性和準(zhǔn)確性。文獻(xiàn)[15]通過求取嵌入維數(shù)與延遲時間，把人流量時間序列重構(gòu)為多維度空間矢量，將有規(guī)律的混沌吸引子恢復(fù)出來后再進(jìn)行預(yù)測。文獻(xiàn)[16]應(yīng)用混沌理論，重構(gòu)GPS衛(wèi)星鐘差序列并進(jìn)行了預(yù)測，取得了較好的短期預(yù)測效果。

因此，本文嘗試將混沌理論與 LSTM相結(jié)合，提出1種新的預(yù)測模型。首先對時間序列進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解（empirical mode decomposition，EMD）并降噪，去除序列包含的白噪聲部分，然后求取時間序列的延遲時間，嵌入維數(shù)以及李雅普諾夫指數(shù)，證明GPS高程站心坐標(biāo)時間序列具有混沌特性，并重構(gòu)序列相空間。最后，將相空間每1維特征向量作為LSTM的時間步輸入網(wǎng)絡(luò)進(jìn)行訓(xùn)練。本文方法將預(yù)測結(jié)果的均方根誤差（root mean square error，RMSE）與平均絕對值百分比（mean absolute percentage error，MAPE）與常用的 LSTM 和 LSTM+卷積神經(jīng)網(wǎng)絡(luò)（convolutional neural networks，CNN）的方法，以及文獻(xiàn)[6，10]提及的方法相比較，希翼通過著重對GPS高程時間序列特性進(jìn)行分析，可彌補(bǔ)簡單應(yīng)用神經(jīng)網(wǎng)絡(luò)可解釋性較差的缺點(diǎn)，對于時間序列變化趨勢的預(yù)測會更為顯著與準(zhǔn)確，以用于GPS高程時間序列的趨勢預(yù)測或補(bǔ)全工作。

1 GPS高程時間序列預(yù)測方法

本文預(yù)測方法的網(wǎng)絡(luò)結(jié)構(gòu)如圖1所示，主要由時間序列降噪、混沌特性分析，以及LSTM預(yù)測網(wǎng)絡(luò)幾部分組成。其中，實(shí)驗(yàn)數(shù)據(jù)來源于國際全球衛(wèi)星導(dǎo)航定位系統(tǒng)服務(wù)組織（The International Global Navigation Satellite Systems（GNSS） Service，IGS）基準(zhǔn)站上海佘山站（SHAO）的2 380 d站心坐標(biāo)系高程數(shù)據(jù)，將其作為原始時間序列。

圖1 GPS高程時間序列預(yù)測模型

1.1 序列降噪

EMD方法是1種數(shù)據(jù)驅(qū)動的降噪方法，認(rèn)為瞬時頻率如果有意義，函數(shù)必須是對稱的，局部均值為零，并且具有相同的過零點(diǎn)和極值點(diǎn)，由此可將原始時間序列分解為若干個包含不同時間尺度的本征模函數(shù)（intrinsic mode function，IMF），如圖2所示。其中包括高頻噪聲項(xiàng)、周期項(xiàng)、以及趨勢項(xiàng)等[17-18]，再通過計(jì)算各IMF與原序列的相關(guān)性重構(gòu)降噪后的時間序列，去除高頻噪聲，各IMF與原始時間序列的相關(guān)系數(shù)如表1所示，實(shí)現(xiàn)對GPS時間序列的降噪后，得到的序列尺寸應(yīng)為2 380×1，即降噪后時間序列共有2 380 d數(shù)據(jù)，每1天數(shù)據(jù)包含1個特征維度。

表1 各本征模函數(shù)與原序列的相關(guān)系數(shù)

圖2 原始時間序列與經(jīng)EMD分解后的各本征模函數(shù)

由圖2可知，原始時間序列經(jīng)EMD分解得到 了10個不同時間尺度的本征模函數(shù)，頻率依次降低，其中IMF10頻率最低且具有明顯的趨勢項(xiàng)特征。而由表1可知第1個相關(guān)系數(shù)極小值對應(yīng)的IMF為IMF4，說明此時IMF序列與原序列相關(guān)性最低，可知IMF1~IMF4代表了高頻噪聲[17-18]；所以將它們歸為噪聲序列，IMF5~IMF10則重構(gòu)為降噪后時間序列，序列降噪前后對比如圖3所示。

由圖3可知，降噪后的時間序列較好地保留了它的周期項(xiàng)特征與趨勢項(xiàng)特征，時間序列的變化情況更為清晰。

1.2 混沌特性分析

降噪后的時間序列仍為1維的時間序列，此時序列的特征并不明顯，需要對序列進(jìn)行合適的采樣，拓展特征維數(shù)，故采用互信息法求取合適的延遲時間，CAO方法求取嵌入維數(shù)，最終得到李雅普諾夫指數(shù)，驗(yàn)證序列具有混沌特性并將降噪后的時間序列重構(gòu)為包含 4個維度特征的多特征時間序列，重構(gòu)后的序列尺寸應(yīng)為2 340×4，即重構(gòu)后的序列包含2 340 d數(shù)據(jù)，每1天的數(shù)據(jù)包含4個特征維度。

1）延遲時間的求取?；バ畔⒎ㄊ?種成熟的確定延遲時間的方法，它能夠描述序列具有的非線性特征，最大化時間序列包含的信息量，其計(jì)算公式可表示為

式中：E(τ)表示時間序列平均信息熵；pi與 pj分別表示時間序列在區(qū)間i、j內(nèi)的概率；pij(τ)為其聯(lián)合概率；τ為延遲時間。由式（1）可知，隨著τ的增大，時間序列 x(t)與 x(t+τ)的 E(τ)會逐漸降低，直到某個臨界點(diǎn)，它們將不再關(guān)聯(lián)，而這個臨界點(diǎn)通常取第1個極小值點(diǎn)。針對降噪后序列應(yīng)用互信息法求得的不同延遲時間下的平均信息熵如圖 4所示。

由圖 4可知，平均信息熵在延遲時間為 14 d的時候，達(dá)到第1個極小值，說明此時序列片段間關(guān)聯(lián)性達(dá)到最低，故選擇延遲時間為14 d。

2）嵌入維數(shù)的求取。CAO方法求取嵌入維數(shù)，能夠使重構(gòu)后的時間序列相空間不存在信息折疊，在混沌時間序列分析時廣泛使用[19]，算法步驟如下：

步驟1：從較小的嵌入維數(shù)m開始，分別重構(gòu)出m維與m+1維的相空間，計(jì)算相空間中鄰近點(diǎn)對間的歐式距離比值()，rim，可表示為

步驟3：計(jì)算E(m)與E(m+1)的變化比Em，可表示為

步驟4：重復(fù)步驟1~步驟3直到Em穩(wěn)定。

圖4 不同延遲時間下的平均信息熵

若序列存在混沌特性，那么當(dāng) m增加到某個值時，Em會收斂于1，表明此時時間序列的吸引子在相空間已完全打開，嵌入維數(shù)的增加不會再使特征序列包含的信息量變多，所以將 m+1作為嵌入維數(shù)最合適，根據(jù)CAO方法求得的各維數(shù)下Em的變化情況如圖5所示。

圖5 最小嵌入維數(shù)

由圖5可知，維數(shù)達(dá)到4以后，變化比在小范圍內(nèi)逐漸穩(wěn)定，可以認(rèn)為在該維度上重構(gòu)時間序列能夠很好地展開11維時間序列，近似反映原系統(tǒng)相空間的變化規(guī)律，因此嵌入維數(shù)選擇為4。

3）GPS時間序列的混沌性判別?；煦缧耘袆e是混沌時間序列相空間重構(gòu)的可行性保證，常用方法是通過求取延遲時間和嵌入維數(shù)來確定最大李雅普諾夫指數(shù)，進(jìn)而判斷時間序列是否具有混沌特性。

李雅普諾夫指數(shù)反映了混沌吸引子中鄰近軌道的分離速度，而最大李雅普諾夫指數(shù)決定了相空間相鄰軌道是否穩(wěn)定，若其大于 0，則序列變化會以指數(shù)增長，直到無法預(yù)測。根據(jù)求得的延遲時間與嵌入維數(shù)計(jì)算李雅普諾夫指數(shù)，具體步驟如下：

步驟 1：對降噪時間序列做快速傅里葉變換，并通過周期-功率譜圖估算序列平均周期；

步驟2：利用延遲時間τ，嵌入維數(shù)m重構(gòu)時間序列，序列x(n)經(jīng)過重構(gòu)后可表示為

式中： xn?m+1表示重構(gòu)后的第 n-m+1維時間序列；x(n ?m+1)表示原始序列第n-m+1個值。

步驟3：找到每個點(diǎn) x (t)的最近鄰點(diǎn)x(t?)，同時限制短暫分離，計(jì)算鄰近點(diǎn)對的距離初值 xt(0)，可表示為

式中：p 為平均周期；t=1，2，…，n-m+1。

步驟 4：計(jì)算每個點(diǎn)χ（t）在 i個延遲后相應(yīng)鄰近點(diǎn)對的距離χt（i），即

式中：t?=1，2，… ，m，且 t≠t?；i=1，2，… ，min(m-t，m-t?)。對每個 i求取所有 t的χt（i），并求其平均χ（i），在χ（i）與i有線性關(guān)系的區(qū)域用最小二乘法擬合直線，該直線斜率即為最大李雅普諾夫指數(shù)。

根據(jù)上述步驟求取時間序列的平均周期，時間序列周期功率譜如圖6所示。

圖6 周期功率譜（即序列長度時）

由圖6可知，序列周期圖為2 380 d時，功率達(dá)到最大點(diǎn)，其次為平均周期取793 d時，故取793 d為序列平均周期，結(jié)合延遲時間與嵌入維數(shù)求取最大李雅普諾夫指數(shù)的結(jié)果如圖7所示。

圖7 時間序列最大李雅普諾夫指數(shù)

根據(jù)圖7可知，實(shí)線為擬合的直線，其斜率為0.033 4，即最大李雅普諾夫指數(shù)大于0，證明實(shí)驗(yàn)所用高程時間序列具有混沌特性，可以采用混沌理論重構(gòu)GPS高程時間序列相空間并進(jìn)行預(yù)測。

1.3 LSTM預(yù)測模型

LSTM是1種特殊的循環(huán)神經(jīng)網(wǎng)絡(luò)，有著與循環(huán)神經(jīng)網(wǎng)絡(luò)相似的單元結(jié)構(gòu)，能夠通過前向與反饋傳播有效記憶時間序列時間節(jié)點(diǎn)間的變化特征。除此之外，LSTM增加了遺忘門、輸入門和輸出門。其中遺忘門能夠限制上一時間節(jié)點(diǎn)輸出的信息量，輸入門可以選擇性地獲取新的輸入信息，輸出門則能根據(jù)新輸入的信息與上一時刻輸出的信息計(jì)算出當(dāng)前時刻的輸出；3個門結(jié)構(gòu)相互作用，通過反復(fù)訓(xùn)練獲得網(wǎng)絡(luò)權(quán)值，它能夠充分學(xué)習(xí)到不同時間步之間短期或長期的時間特征。

為加快網(wǎng)絡(luò)權(quán)值的收斂，故先對重構(gòu)的時間序列進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化處理，再輸入 LSTM預(yù)測模型進(jìn)行學(xué)習(xí)，并計(jì)算網(wǎng)絡(luò)輸出的迭代預(yù)測值與真值的損失（loss），使用 ADAM（adaptive moment estimation）方法對網(wǎng)絡(luò)權(quán)重值不斷優(yōu)化改正，最終得到預(yù)測模型。

1.4 評價指標(biāo)

預(yù)測結(jié)果的評價指標(biāo)采用 RMSE和 MAPE，其中：RMSE對異常點(diǎn)較敏感，本文用它刻畫真實(shí)值與預(yù)測值的偏差，當(dāng)預(yù)測時出現(xiàn)的不合理點(diǎn)越少，RMSE也會越?。籑APE則考慮了誤差與真實(shí)值之間的比例，用它來反映預(yù)測值與真實(shí)值的相對誤差，它們的計(jì)算公式可分別表示為

式中：n 為序列長度；true［i］與 predict［i］分別代表序列第i個真實(shí)值和預(yù)測值。

本文還將預(yù)測點(diǎn)與真實(shí)值的趨勢擬合情況作為評價指標(biāo)之一，并且分析了本文方法應(yīng)用在其他站點(diǎn)數(shù)據(jù)上的泛化能力。

2 實(shí)驗(yàn)與結(jié)果分析

本文的實(shí)驗(yàn)數(shù)據(jù)來源于中國地震局全球衛(wèi)星導(dǎo)航定位系統(tǒng)（global navigation satellite systems，GNSS）數(shù)據(jù)產(chǎn)品服務(wù)平臺，共設(shè)計(jì)了2組實(shí)驗(yàn)用于驗(yàn)證本文提出的預(yù)測模型的有效性。第1組實(shí)驗(yàn)基于IGS基準(zhǔn)站上海佘山站（SHAO）的2 380 d站心坐標(biāo)系高程數(shù)據(jù)，對比方法為目前常用的LSTM和 LSTM+CNN方法、文獻(xiàn)[6]和文獻(xiàn)[10]使用的變分模態(tài)分解（variational mode decomposition，VMD）以及廣義回歸神經(jīng)網(wǎng)絡(luò)（generalized regression neural network，GRNN）的預(yù)測方法。第 2 組實(shí)驗(yàn)針對SHAO、武漢站（WUHN）、拉薩站（LHAS）、新疆鄯善站（XJSS）和北京房山站（BJFS）等5個不同區(qū)域站點(diǎn)的各2 380 d站心坐標(biāo)系高程數(shù)據(jù)，分別使用不同方法進(jìn)行預(yù)測實(shí)驗(yàn)，旨在驗(yàn)證本文預(yù)測模型的泛化能力。實(shí)驗(yàn)中數(shù)據(jù)預(yù)處理以及混沌特性分析部分的實(shí)驗(yàn)是基于MATLAB2014完成的，預(yù)測模型的搭建與實(shí)驗(yàn)是基于 Python3.7與Tensorflow1.9完成的。

2.1 預(yù)測模型的訓(xùn)練

綜合前文的分析，結(jié)合圖1可知，重構(gòu)后的時間序列共有2 340 d，每1天的數(shù)據(jù)有4個維度，實(shí)驗(yàn)選取了前2 310 d數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，剩余部分作為測試數(shù)據(jù)，并且通過網(wǎng)格搜索法選取預(yù)測模型超參數(shù)，最終選定隱層數(shù)量為32層，每1隱層中的LSTM神經(jīng)元為30個，訓(xùn)練輪數(shù)為500輪，學(xué)習(xí)率為0.001，批尺寸為64個。模型的輸入為重構(gòu)后序列的第1~2 309 d數(shù)據(jù)，與之對應(yīng)的輸出標(biāo)簽為原始時間序列中第2~2 310 d數(shù)據(jù)，訓(xùn)練過程中各方法的損失函數(shù)收斂情況如圖8所示。

由圖8可知，在訓(xùn)練集上，所有方法的loss收斂都很快。其中LSTM方法的loss最大，收斂后有一定的波動，LSTM+CNN的loss較之更低，且較穩(wěn)定。本文方法與文獻(xiàn)[10]方法的 loss大小明顯小于前 2種方法，說明在訓(xùn)練集上更容易學(xué)到序列的變化模式，效果優(yōu)于前2種方法；不同在于本文方法loss波動較大，但總體看來 loss略低于文獻(xiàn)[10]的方法，說明本文方法在訓(xùn)練集上表現(xiàn)最優(yōu)，同時也印證了模型超參數(shù)的選擇是較為合理的。

圖8 損失函數(shù)收斂

2.2 SHAO站預(yù)測結(jié)果對比

考慮到迭代預(yù)測的結(jié)果存在一定的不確定性，所以本文基于4種方法分別訓(xùn)練了20次模型，并分別統(tǒng)計(jì)了預(yù)測結(jié)果的評價指標(biāo)，見圖9、圖10，各方法最優(yōu)結(jié)果的評價指標(biāo)如表2所示。

表2 各方法評價指標(biāo)對比

綜合圖9、圖10與表2可知，所有方法的預(yù)測結(jié)果都在一定的范圍內(nèi)浮動，是相對穩(wěn)定的。其中，文獻(xiàn)[6]由于通過 Lag Feature方式重構(gòu)時間序列，并進(jìn)行了網(wǎng)絡(luò)超參數(shù)優(yōu)選，因此 RMSE相較前2種方法有明顯提高。文獻(xiàn)[10]則結(jié)合了混沌理論重構(gòu)的多特征時間序列和 GRNN方法，同時提高了RMSE和MAPE，說明了多特征預(yù)測的有效性。本文方法同樣利用混沌理論構(gòu)建多特征，但采用了 LSTM進(jìn)行預(yù)測：得到的 MAPE為2.82 %，說明預(yù)測值與真值相對誤差?。欢鳵MSE僅為0.001 2，說明預(yù)測結(jié)果中異常點(diǎn)較少，其對應(yīng)的觀測值物理意義為站心坐標(biāo)預(yù)測值與真值間的平均誤差為 1.2 mm，這能夠滿足 GPS靜態(tài)數(shù)據(jù)后處理的精度要求、在所列方法中，預(yù)測效果為最佳。各方法對 SHAO站 30 d最佳的滾動預(yù)測效果圖見圖11。

由圖 11可知：LSTM方法和 LSTM+CNN方法預(yù)測點(diǎn)擬合的整體趨勢較為平緩，對于序列的局部變化難以及時做出調(diào)整進(jìn)行預(yù)測；而文獻(xiàn)[10]和本文方法由于使用了多特征，使得預(yù)測結(jié)果對于單調(diào)性以及凹凸性發(fā)生變化的點(diǎn)更為敏感，預(yù)測的趨勢更為準(zhǔn)確。這也說明了為什么本文方法在訓(xùn)練時，雖然誤差總體較小，但波動較大。

2.3 模型的泛化能力分析

為說明本文提出的預(yù)測方法具有較好的泛化能力，本文針對WUHN站、LHAS站、BJFS站、XJSS站分別做了4種方法的預(yù)測實(shí)驗(yàn)，預(yù)測效果如圖12所示。

圖9 均方根誤差（RMSE）對比

圖10 平均絕對值百分比（MAPE）對比圖

圖 11 SHAO 站 30 d 預(yù)測效果

圖12 各站30 d的預(yù)測效果圖

由圖12可知，針對不同地區(qū)的高程觀測數(shù)據(jù)，各站點(diǎn)的預(yù)測結(jié)果與 SHAO站大致相同，LSTM、LSTM+CNN的預(yù)測測結(jié)果較為平緩，能夠反映一定的整體趨勢，而本文的方法不僅能準(zhǔn)確反映序列整體的變化趨勢，局部點(diǎn)的變化趨勢也更為顯著，說明該預(yù)測模型具有一定的泛化能力，同時也說明了混沌分析，用多特征變量進(jìn)行預(yù)測以及建立混合模型的有效性。

基于上述實(shí)驗(yàn)結(jié)果，可以得出以下結(jié)論：

1）單變量特征的預(yù)測方法在精度和預(yù)測曲線形態(tài)的擬合上要略差于多變量特征預(yù)測方法；

2）基于混沌理論將單一維度時間序列重構(gòu)為多特征時間序列，是1種有效的拓展或恢復(fù)特征的方法，能夠有效提高預(yù)測的效果；

3）結(jié)合對時間序列自身特性分析與預(yù)測模型的混合模型預(yù)測方法要優(yōu)于某一預(yù)測網(wǎng)絡(luò)模型的簡單應(yīng)用。

3 結(jié)束語

針對GPS高程時間序列的特點(diǎn)，將混沌理論與LSTM相結(jié)合，本文提出了1種新的混合預(yù)測模型。首先對時間序列進(jìn)行EMD降噪，去除白噪聲等無關(guān)信息；然后驗(yàn)證了GPS高程時間序列具有混沌特性；并將1維時間序列在高維展開，增加了特征維度，最大化序列包含信息；再作為LSTM的時間步輸入模型進(jìn)行訓(xùn)練。通過 2組實(shí)驗(yàn)數(shù)據(jù)計(jì)算結(jié)果表明，本文方法的預(yù)測結(jié)果中，RMSE最低可達(dá) 0.001 2，而 MAPE僅為2.9 %，能夠滿足GPS后處理的精度要求，提高了對高程時間序列的預(yù)測精度以及變化趨勢的準(zhǔn)確性，同時具有一定的泛化能力，是深度學(xué)習(xí)在GPS數(shù)據(jù)分析領(lǐng)域的 1次應(yīng)用，可用于今后的實(shí)際工程中的預(yù)測或補(bǔ)全應(yīng)用。接下來的研究方向會針對訓(xùn)練時不穩(wěn)定的現(xiàn)象進(jìn)行模型改進(jìn)，同時探索深度學(xué)習(xí)在GPS領(lǐng)域的其他應(yīng)用。