BDS衛星鐘差短期預報的LSTM算法

吉長東，朱錦帥,，黎 虎，張 萌，呂廣涵

（1. 遼寧工程技術大學 測繪與地理科學學院，遼寧 阜新 123000；2. 張家界市測繪院，湖南 張家界 427000）

0 引言

在多導航系統并存的環境下，實時衛星鐘差產品的精度直接影響到全球衛星導航系統（global navigation satellite system, GNSS）的服務能力[1-2]。我國北斗衛星導航系統（BeiDou navigation satellite system, BDS）搭載的原子鐘均為自主研發，其精度和穩定性與美國全球定位系統（global positioning system, GPS）和歐盟伽利略衛星導航系統（Galileo navigation satellite system, Galileo）的衛星鐘不同，因此，BDS鐘差產品的精度和穩定性應受到更多的關注[3]。

超快速鐘差產品作為國際GNSS服務組織（International GNSS Service, IGS）的核心產品之一，是實現分米級實時定位服務的重要一環。文獻[4-5]通過二次多項式加1個周期項擬合，預報IGS超快速產品中的GPS鐘差，但其精度與廣播星歷相比并沒有顯著改善。文獻[6]把小波神經網絡（wavelet neural network, WNN）模型應用到鐘差短期預報中，將小波分析與神經網絡結合，預報精度優于二項式模型和灰色模型。文獻[7]以WNN模型為基礎，對數據預處理方法進行改進后，以IGS發布的GPS超快速鐘差進行實驗并得到優于超快速鐘差預報部分數據的結果。文獻[8]使用徑向基函數（radial basis function, RBF）神經網絡，對GPS鐘差進行超短期5 min和1 h以及短期1 d預報，均可取得高精度預報結果。文獻[9]提出的循環神經網絡（recurrent neural network, RNN）是1種具有記憶功能的深度學習模型。文獻[10]通過修改隱藏層的結構，建立1種基于RNN的長短時記憶（long short-term memory, LSTM）模型，進而削弱了RNN在長序列訓練過程中的梯度消失和梯度爆炸問題。

目前，只有德國地學研究中心（Deutsches Geo Forschungs Zentrum, GFZ）和我國的全球連續監測評估系統（international GNSS monitoring and assessment system, iGMAS）提供開放式BDS鐘差產品，且現有的鐘差預報算法主要集中于穩定的GPS鐘差數據以及IGS發布的鐘差產品，對于BDS以及iGMAS發布的鐘差產品缺少研究。為此，針對BDS超快速鐘差，提出1種將修正1次差的鐘差預處理方法和LSTM模型相組合的鐘差短期預報方法。

1 鐘差數據預處理及LSTM預報原理

iGMAS綜合鐘差產品的時標采用北斗時（BDS time, BDT），用戶將iGMAS的BDT產品轉換到GPS時（GPS time, GPST）時，不考慮微小量可以時[11]，其計算方法為

式中：涉及物理量的單位均為 s；tGPST為t時刻GPS鐘差值；tBDT為t時刻BDS鐘差值。iGMAS提供的鐘差產品具體數據如表1所示，鐘差產品精度用均方根值（root mean square，RMS）表示。

表1 iGMAS鐘差產品介紹[11]

1.1 鐘差1次差預報方法原理

為了獲得高精確度的時間信息值，原始衛星鐘差相位數據具有更多的有效位數。文獻[6]提出了對鐘差相位數據相鄰歷元間求1次差的方法，這不僅可以降低原始鐘差數據序列中趨勢項的影響，進而得到1組有效數字位數減少且數值較大的1組數據序列。設L=｛li,i=1,2,…,n｝（li為i時刻鐘差值）為任意1組鐘差序列，對序列中2個相鄰歷元的鐘差做差，可以得到鐘差的1次差分序列ΔL=｛Δli,i=2,3,…,n｝，其中Δli=li-li-1。利用1次差分后的序列ΔL進行學習、建模并預報第n個歷元以后m（m＞0）個歷元的鐘差的1次差分序列ΔLm=｛Δlj,j=n+1,n+2,…,n+m｝，最后根據預報的1次差序列求累加和與第n個歷元的鐘差值ln相加，即可得第k個歷元的鐘差相位值[6]為

式中：Δlk為預報的第k個歷元的1次差值，k=n+1,n+2,…,n+m。

1.2 基于鐘差1次差分序列的粗差處理

針對鐘差相位數據相鄰歷元間求1次差后的數據序列，文獻[12]設計了1種修正1次差法數據預處理方法，以改進中位數為基礎的抗差估計探測粗差方法：將每個鐘差1次差分數據Δli與1次差分序列的中數及絕對中位數（median absolute deviation，MAD）數倍之和進行比較，MAD可表示為

式中：kΔ為1次差序列的中間數；Median為對序列取中位數,即kΔ=Median｛ΔL｝。若鐘差1次差分值滿足式（4），則將其判定為異常值并剔除該值，然后通過內插或者置零的方法補充該點數據[13]，即

式中：a為根據實際情況需要確定的1個正整數。

本文遵循探測出的異常值個數不超過建模數據的5%的原則，將a設置為3，這樣既在一定程度上降低異常值對模型預報性能的影響，又避免了有效信息的剔除。對于剔除數據歷元，這里使用線性插值的方法補全歷元。

1.3 循環神經網絡原理

循環神經網絡是1種可以有效提取、利用和處理時間序列的高位非線性動力學模型[14]。在循環神經網絡中，神經元不但可以接受其他神經元的信息，也可以接受自身的信息，將神經元最近1次（或幾次）的結果導入自己接下來的運算中，形成具有環路的網絡結構。設x=(x1,x2,xt)為模型輸入，y=(y1,y2,yt)為模型輸出，在1個完全連接的RNN模型中，在時刻t時，令xt表示網絡的輸入，yt表示網絡輸出，網絡模型如圖1所示。

圖1 RNN模型及其展開示意圖

圖1中：h為隱藏層狀態（即隱藏層神經元的活性值，也稱之狀態或隱狀態），ht-1、ht-1為t-1、t時刻活性值；U、W和V均為網絡參數，U為狀態—狀態的權重矩陣，W為狀態—輸入的權重矩陣，V為狀態—輸出的權重矩陣。則RNN模型在時刻t的參數更新公式為：

式中：zt為t時刻隱藏層的凈輸入；b表示偏置向量；f(zt)表示非線性的激活函數。從圖1和式（5）、式（6）可以看出，t時刻隱藏層的狀態ht不僅有與時刻t的輸入xt相關，也和上一時刻隱藏層的狀態ht-1有關。

1.4 LSTM原理

在RNN實際應用中，隱藏層存儲的時間是有限的，同時梯度消失的問題也是長程時間序列預測中存在的主要問題。基于此，文獻[15]以LSTM模型為基礎，在網絡中引入1個新的內部狀態ct和門控機制，專門進行線性的循環信息傳遞，同時（非線性的）輸出信息給隱藏層的外部狀態th。在t時刻，內部狀態tc的公式計算為：

式中fsig為sigmoid型激活函數。

2 鐘差預報的LSTM模型構造

利用LSTM模型建立鐘差預報模型的方式包括鐘差數據的訓練和預報兩部分。在訓練階段，將鐘差序列作為神經網絡的訓練數據，輸入的訓練數據既作為網絡的輸入也作為網絡的期望輸出，通過LSTM神經元的訓練，使其掌握數據前后間的映射關系；在預報階段，將當前時刻的鐘差數據作為網絡的輸入，通過已建立的網絡之間前后映射及拓撲關系，依次進行外推得到預報結果。

根據訓練輸入的鐘差序列有限個樣本點的數據特征，參考LSTM預報模型的從簡設計原則，構建基于LSTM的深度循環神經網絡。對于鐘差序列預報的整體框架如圖2所示。

LSTM模型預報部分由輸入層、隱藏層、輸出層、網絡訓練和參數優化5個部分組成。輸入層負責對輸入的鐘差時間序列進行處理，使得輸入的鐘差時間序列滿足LSTM神經元的輸入要求；隱藏層使用基于LSTM神經元的深度循環神經網絡；網絡訓練部分使用自適應優化算法中的亞當（Adam）算法結合LSTM神經元對輸入數據進行訓練；使用網格搜索和交叉驗證算法進行超參數的優化；輸出層負責預報結果的輸出。

圖2 LSTM預報模型流程

以預處理后鐘差時間序列L0= {l1,l2,…,lr}為例，將其導入網絡進行訓練、預報計算的具體流程如下：

1）數據標準化。將預處理后的鐘差時間序列L0= {l1,l2,…,lr}分 為 訓 練集Ltr= {l1,l2,…,ls}和 測試集Lte= {l s+1,ls+2,…,lr}，其中1＜s＜r。之后，對訓練集中的每個數據lu(1 ≤u≤s)通過式（14）進行Z-score標準化，得到標注化數據集Lt′r= {l1′,l2′ ,…,ls′}，即

式中：μ為序列Ltr的均值；σ為序列Ltr的標準差。

2）數據分割。根據數據分割長度L，將標準化后的訓練數據集Lt′r分割為長度為L的(s-L+1)個數列，使訓練集輸入變為X= {X1,X2,…,Xs-L+1}，其 中Xp= {l′p,l′p+1,…,l′p+L-1}；則此刻 對應 的 理論 輸出Y= {Y1,Y2,…,Ys-L+1}，其中Yp={l′p+1,l′p+1+2,…,l′p+L}，1≤p≤s-L+1。

3）數據訓練。當模型的輸入為X時，令其經隱 藏 層 計 算 后 的 輸 出 為P= {P1,P2,…,Ps-L+1}，則Pp=fLSTM-forward(X p,c p-1,hp-1)，Xp表示當前時 刻 的輸入，cp-1表示前一歷元LSTM神經元的狀態值，hp-1表示前一歷元LSTM神經元的輸出值，fLSTM-forward表示LSTM神經元的更新公式。網絡訓練過程中以求得損失函數結果更小為模型優化目標，在這個過程中使用小批量梯度下降法（minibatch gradient descent, M-BGD），每次選擇批尺寸大小為B個輸入值作為1批輸入網絡進行訓練，通過Adam優化器使得網絡不斷更新各個參數的權重，最后得到1個對輸入序列X擬合效果最好的LSTM模型，記作fLSTM-net。

4）超參數優選。將網格搜索和交叉驗證方法結合，對數據分割長度L和批尺寸大小B進行優選。

5）數據預測。預測是采用逐個歷元數據迭代的方法，利用訓練好的fLSTM-net模型，將訓練集的最后1組 數 據輸 入fLSTM-net中，得 到 的 輸 出 結 果 為Ps-L+1=fLSTM-net(Ys-L+1)={ps-L+2,ps-L+3,…,ps+1}。將s+1時刻的預測值ps+1與YL合 并 為1組新 的 數據將YL+1帶入fLSTM-net中，可以得到s+2時刻的預測值ps+2；按照這個方法依次計算，得到預測序列P0={p s+1,ps+2,…,pr}。將預測序列P0進行Z-score反標準化還原，就可以得到預報結果。

3 實驗與結果分析

實驗數據選取iGMAS超快速（iGMAS ultrarapid, ISU）鐘差文件，其中包括超快速觀測部分（ISU observed, ISU-O）和超快速預報部分（ISU predicted, ISU-P），數據文件格式為isu*****_00.clk，*****為BDT，前4位是BDT星期，第5位是星期內天）。考慮BDS包含地球靜止軌道（geostationary Earth orbit, GEO）、傾斜地球同步軌道（inclined geosynchronous orbits, IGSO）及中圓地球軌道（medium Earth orbit, MEO）三種不同軌道的衛星，不同軌道所處空間環境不同，因此在每個軌道隨機選取某一刻衛星進行實驗，本文選取的是C04、C07和C12三顆衛星進行實驗。最后，以實測部分ISU-O為訓練數據和期望輸出，采用均方根誤差（root mean squared error，RMSE）和極差（Range）2個性能指標來評估預報結果的精度與穩定性[9]。RMSE和R的具體定義為：

式中：R為極差；εi為i時刻期望輸出li與預報值l?i的殘差；εmax和εmin分別代表差值序列中的最大值與最小值。

3.1 超參數優選實驗分析

結合網格搜索和交叉驗證的方法，對LSTM神經網絡中使用的超參數窗口分割長度L和批尺寸大小B兩種超參數進行優化。根據數據類型設計超參數的取值范圍，其中窗口分割長度分別取0.5、1.0、1.5、2.0、3.0、4.0、6.0及12 h內的數據個數，即L∈{2, 4, 6, …, 24, 48}，批尺寸大小B∈{2, 4,… ,48}，以訓練集中訓練數據的誤差函數最小為目標進行訓練，運用多層網格搜索對超參數L、B進行優選。

以C12衛星的ISU鐘差預報為例，選取2019-01-14的15 min采樣率共96個歷元的ISU-O數據進行訓練，預報96個歷元數據。以2019-01-15的ISU-O數據為真值，計算RMSE值，遍歷本文所選的L、B兩個超參數的取值范圍，繪制超參數優選圖如圖3所示。

圖3 參數優選

從圖3中可以看出，當批尺寸大小B的取值在20以內時，LSTM網絡的預報結果不甚理想，且容易陷入局部最優解，B的取值大于20時，不同的窗口分割長度預報結果的RMSE值變化在0～1.2×10-8的范圍內。表2給出了C04、C07和C12三顆衛星使用上述預報方案時，前5組最優參數組合以及預測精度。綜合來看，當L=4、B=30時，預報精度更為理想。

表2 C04、C07、C12衛星LSTM模型前5組最優參數組合及預測精度

3.2 不同數據預處理方法的LSTM預報實驗

不同數據預處理方法對神經網絡的預報精度有著重要影響，使用文中提到的2種1次差數據預處理方法與未使用1次差處理的數據進行LSTM模型預報精度分析。為了便于區分，將未使用1次差處理鐘差數據的LSTM預報方案命名為LSTM_1，使用鐘差1次差法處理的預報方案命名為LSTM_2，經過改進MAD法粗差處理后的1次差數據預報方案命名為LSTM_3。為了對比這3種方案的預報性能，以日期2019-01-09的ISU-O數據進行預處理，然后訓練預報第2天（2019-01-10）的鐘差，以第2天的ISU-O數據為真值與預報結果作對比，其部分預報結果的相位數據及1 d的預報殘差如圖4所示。

將圖4中實驗結果與對應的真值進行對比后可以看出，LSTM_2方案的預報殘差變化范圍最大，其前5個歷元的預報結果最好，經過1次差處理后，鐘差數據中的粗差被放大，導入網絡訓練的輸入數據中粗差對網絡訓練的影響增大，干擾了模型的訓練、擬合過程，造成網絡出現過擬合現象，使網絡結構極不穩定，因此當預報數據長度超過20個歷元后，其預報殘差將持續增大（包括向正無窮與負無窮）。對比LSTM_1和LSTM_3這2種實驗方案：前5個歷元，2種方法的殘差大小和變化趨勢有明顯差距，在前20個歷元內，2種方法的預報殘差趨于穩定，變化趨勢相近；但是LSTM_3的殘差波動幅度明顯較LSTM_1的殘差浮動小。

圖4 不同預處理方法后預報結果與殘差圖

為進一步對比3種預處理方法對預報結果精度的影響，將預報結果與ISU-O真值進行RMSE和R值統計，如表3所示。

從表3可以看出：LSTM_2由于原始數據中粗差被放大造成網絡訓練過程中不穩定，其R值從6 ns到30 ns，而RMSE值高達16 ns，預報結果的精度和穩定度較差；對比LSTM_1和LSTM_3，這2個模型的RMSE和R值均在納秒級，其中以LSTM_3方案的結果最好，其RMSE和R值都能保持在2 ns以內。通過對1次差序列進行改進MAD法粗差探測、剔除與補齊，不僅削弱了鐘差相位數據中趨勢項和粗差數據對網絡訓練的影響，還能減少輸入數據有效數位長度，使輸入的數據更能表征數據的變化特點，減少網絡的冗余。

表3 不同數據預處理方法預報精度統計表 單位：ns

3.3 ISU鐘差預報精度分析

圖5 ISU-O短期預報結果與殘差

為了進一步驗證LSTM在BDS的ISU鐘差產品中的預報性能，選用日期2019-01-09的ISU-O數據作為LSTM和WNN模型的訓練數據，以ISU鐘差產品中的預報部分和WNN預報結果作為對照。選擇3.2實驗中LSTM_3使用的數據預處理方法對實驗數據進行預處理，然后將預處理后的數據，分別通過WNN和LSTM兩種神經網絡模型進行訓練并預報第2天（2019-01-10）的鐘差，以第2天的ISU-O數據為真值，對比ISU-P數據以及WNN和LSTM預測結果，其部分預報結果的相位數據及1 d的預報殘差如圖5所示。

從圖5每顆衛星前12個歷元的預報鐘差值可以發現：LSTM、WNN和ISU-P的結果與ISU-O的數據趨勢基本相同，結合預報殘差圖分析，ISU-P數據的殘差偏離程度最大，尤其是C04衛星，其ISU-P數據的殘差呈遞增趨勢，其預報誤差達到了30 ns以上；WNN的預報結果具有明顯的趨勢項，從C04和C07 2顆衛星的預報結果來看，其預報誤差呈明顯的遞增或遞減趨勢，但是其離散程度與ISU-P的殘差相比明顯較小；而LSTM的預報結果較前二者明顯更好一些，表現了明顯的“自我約束能力”。

為了反映LSTM、WNN和ISU-P短期預報的性能，將3種預報結果前1、2、3、4、6及12 h的預報結果精度進行RMSE統計，如表4所示。

表4 LSTM、WNN和ISU-P短期預報RMSE統計表

從表4可以看出：3種預報方案在最初1～2 h的預報精度波動較為明顯，隨著預報時間增加，預報精度逐步提升。為了整體了解LSTM、WNN和ISU-P 3種結果的預報精度和穩定度，將3顆衛星的3種結果與ISU-O作比較表征其精度與穩定度的RMSE與R值進行統計如表5所示。

表5 LSTM、WNN和ISU-P 1d預報精度統計表 單位：ns

從表5中可以看出，LSTM模型的預報精度和穩定性在3顆衛星中表現得最好，3顆衛星平均精度在1 ns左右，穩定度在2 ns左右，尤其是C12衛星的預報精度可達0.85 ns，穩定度在1 ns左右。而ISU-P的預報精度和穩定性最差，1 d內誤差最大可達33 ns，同時精度也達到18 ns，這樣的精度在進行單點定位時誤差在米級。

對比這3種預報結果，WNN和LSTM的日預報精度均優于ISU-P數據，其中在C04、C07和C12三顆衛星中，LSTM的預報精度較ISU-P數據相比精度分別提高了93%、59%和70%，較WNN分別提高了45%、50%和38%。相較于WNN模型，LSTM更適合作為1種較有效的鐘差預報方法應用到鐘差短期預報中。

4 結束語

為了提高鐘差短期預報性能，提出1種1次差和LSTM模型組合的鐘差預報模型。經過超參數優選及2種數據預處理方法對比，搭建適合BDS超快速鐘差預報的組合LSTM模型，將LSTM預報結果分別與WNN模型預報結果及ISU-P數據比較，得到如下結論：

1）采用網格搜索和交叉驗證結合的方法對L和B共2個超參數進行優選，提高了LSTM模型的預報性能和泛化能力。當L=4、B=30時，LSTM模型預報精度較好。

2）對鐘差1次差分數據進行粗差剔除，不僅可以消除鐘差數據中的異常點，還減少了有效數據長度，可降低LSTM網絡的復雜程度。通過這種預處理方法后，LSTM模型預報結果較直接使用LSTM模型進行預報的精度提高了76%。

3）鐘差數據1次差剔除粗差預處理后，再使用LSTM模型預報鐘差，其單日預報精度較WNN預報模型提高了45%，與ISU-P數據相比，則提高了74%。該模型至少在日內預報精度優于ISU-P數據，可以作為1種較為有效的鐘差預報模型。