楊 旭，何祥祥，程茂原，鐘東升，葛雅倩

(1.安徽理工大學(xué) 礦區(qū)環(huán)境與災(zāi)害協(xié)同監(jiān)測煤炭行業(yè)工程研究中心，安徽 淮南 232001；2.安徽理工大學(xué) 空間信息與測繪工程學(xué)院，安徽 淮南 232001；3.安徽理工大學(xué) 礦山采動災(zāi)害空天地協(xié)同監(jiān)測與預(yù)警安徽普通高校重點實驗室,安徽 淮南 232001；4.河海大學(xué) 地球科學(xué)與工程學(xué)院，南京 211100)

在衛(wèi)星導(dǎo)航系統(tǒng)中，時間測量是位置測量的基礎(chǔ)，高精度的星載原子鐘預(yù)報起著重要的作用，將會影響高精度PNT(定位,導(dǎo)航和授時)服務(wù)的質(zhì)量。目前我國的北斗三號衛(wèi)星導(dǎo)航系統(tǒng)已經(jīng)全面建成，開展其高精度預(yù)報有助于提升我國北斗系統(tǒng)的創(chuàng)新服務(wù)能力。

為了滿足實時和近實時用戶對于鐘差時效性的要求，國內(nèi)外專家提出了多種類型的鐘差短期預(yù)報(1 d以內(nèi))方法，主要集中在線性、二次多項式(QP)、譜分析、灰色(GM)、神經(jīng)網(wǎng)絡(luò)(BP)、灰色神經(jīng)網(wǎng)絡(luò)、小波神經(jīng)網(wǎng)絡(luò)(WNN)、自回歸、卡爾曼濾波以及長短時記憶網(wǎng)絡(luò)(LSTM)等模型[1]。不同的預(yù)報模型，其適用范圍有所不同。QP模型添加了衛(wèi)星鐘的物理特性，因而在廣播星歷預(yù)報中廣泛使用，但在中長期預(yù)報中，存在平方項誤差快速增大的問題。GM模型只需要4個起始數(shù)據(jù)即可預(yù)報，崔先強和焦文海提出使用GM模型進行鐘差預(yù)報，實驗結(jié)果表明，GM模型短期預(yù)報與QP模型精度相當，長期預(yù)報中則較QP模型更好[2]。傳統(tǒng)BP模型鐘差預(yù)報中存在陷入局部最優(yōu)和過度學(xué)習的情況，預(yù)報精度時好時壞。陳曉陽基于鐘差一次差分序列，利用LSTM模型進行了衛(wèi)星鐘差短期和中長期預(yù)報，實驗表明，短期預(yù)報中，其預(yù)報性能與QP,GM,BP,WNN相近，中長期預(yù)報中，預(yù)報性能優(yōu)于其他4種模型[3]；S.Wold 和 C.Albano 等人于1983年首次提出偏最小二乘算法(PLS)，它集主成分分析、典型相關(guān)分析和多元線性回歸分析等分析方法的優(yōu)點于一身，是一種新型的多元統(tǒng)計數(shù)據(jù)分析方法，可以有效地克服目前回歸建模中無法解決的樣本容量小于變量個數(shù)時的回歸建模、多個因變量對多個自變量的同時回歸分析等實際問題[4-5]，已在化學(xué)、醫(yī)學(xué)、生物、機械、經(jīng)濟學(xué)等領(lǐng)域中應(yīng)用，但還未在衛(wèi)星鐘差預(yù)報中應(yīng)用。

基于機器學(xué)習算法的衛(wèi)星鐘差預(yù)報不需要考慮鐘差數(shù)據(jù)物理特征，能較好歸納數(shù)據(jù)規(guī)律，直接建立預(yù)報模型就能實現(xiàn)衛(wèi)星鐘差預(yù)報，本文將基于PLS算法進一步拓展機器學(xué)習算法在衛(wèi)星鐘差短期預(yù)報中的應(yīng)用，利用基于鐘差相位數(shù)據(jù)的PLS、基于鐘差相位數(shù)據(jù)一次差分的偏最小二乘(DPLS)、LSTM模型進行BDS-2/BDS-3衛(wèi)星鐘差短期預(yù)報，分析該3種算法在衛(wèi)星鐘差短期預(yù)報應(yīng)用中訓(xùn)練數(shù)據(jù)量最佳確定、算法有效性、精度水平等方面的問題。

1 衛(wèi)星鐘差短期預(yù)報原理

1.1 PLS(partial least square)

PLS是主成分回歸分析(Principal Component Regression,PCR)的改進。PCR未考慮對輸出變量的解釋能力[6]，而PLS在PCR的條件上增加了一個新條件：輸入與輸出應(yīng)具有盡可能大的相關(guān)程度。其數(shù)學(xué)原理為：已知原始數(shù)據(jù)中有m個自變量、n個因變量，則原始數(shù)據(jù)的自變量和因變量分別可以表示為X=(xij)m×N,Y=(yij)n×N(N為樣本點個數(shù)，xij,yij分別為自變量組和因變量組的N次標準化觀測數(shù)據(jù)矩陣元素)。為實現(xiàn)PLS回歸的基本思想，要求u1和v1的協(xié)方差最大，u1、v1分別是矩陣X和Y的第一個主成分，即求解下面的優(yōu)化問題[7]：

max{Cov(u1,v1)}=max[X0w1,Y0c1],

(1)

利用拉格朗日乘數(shù)法求出w1和c1滿足：

(2)

下面介紹PLS算法的建模基本步驟[8]：

1)對原始數(shù)據(jù)進行標準化，提取自變量和因變量的第一對主成分u1和v1。自變量X和因變量Y標準化后得到X0和Y0進而求得主成分對：

(3)

2)建立自變量X0和因變量Y0對u1的回歸方程：

(4)

3)用X1和Y1取代X0和Y0重復(fù)以上步驟，直至所提取的主成分貢獻率超過90%或貢獻率變化低于5%，然后使用取得的主成分進行下一步運算。

4)求關(guān)于主成分的回歸方程。設(shè)最終提出的r個主成分u1，u2，…,ur，滿足：

(5)

將式uk=X0wk(k=1,2，…，r)帶入上式即可得到因變量Y0的PLS回歸方程：

(6)

式中:(cij)m×n為PLS回歸方程的系數(shù)。

5)求原始數(shù)據(jù)的PLS回歸方程。將式(6)中的自變量和因變量標準化還原，得到原始數(shù)據(jù)的PLS回歸方程：

yj=d0j+d1jx1+…+dmjxm,j=1,2，…,n.

(7)

偏最小二乘回歸的成分、殘差矩陣有許多優(yōu)良的性質(zhì)，其中之一是成分之間是相互正交的，這在一定程度上消除了多重線性相關(guān)性。偏最小二乘回歸算法的實質(zhì)是按照協(xié)方差極大化準則，在分解自變量變量數(shù)據(jù)矩陣X的同時，也在分解因變量數(shù)據(jù)矩陣Y，并且建立相互對應(yīng)的解釋隱變量與反映隱變量之間的回歸關(guān)系方程，充分體現(xiàn)了偏最小二乘回歸的基本思想。

1.2 DPLS

通常衛(wèi)星定軌得到的衛(wèi)星鐘差相位數(shù)據(jù)相鄰歷元的數(shù)值變化不大，對PLS模型中相鄰鐘差相位數(shù)據(jù)的差值建立預(yù)報模型(DPLS)，可以降低建立PLS模型數(shù)據(jù)的復(fù)雜性和原始鐘差數(shù)據(jù)序列中趨勢項的影響。設(shè)L={li,i=1,2,…,n}(li為i時刻鐘差值)為任意一組鐘差序列，對序列中兩個相鄰歷元的鐘差做差，可以得到鐘差的一次差分序列ΔL={Δli,i=1,2,…,n}，其中Δli=li-li-1。利用一次差分后的序列ΔL進行學(xué)習、建模并預(yù)報第n個歷元以后m(m>0)個歷元的鐘差的一次差分序列ΔLm={Δlj,j=n+1,n+2,…,n+m}，最后根據(jù)預(yù)報的一次差序列求累加和與第n個歷元的鐘差值ln相加，即可得第k個歷元的鐘差相位值為：

(8)

式中:Δlk為預(yù)報的第k個歷元的一次差值，k=n+1,n+2,…,n+m。

1.3 LSTM

LSTM神經(jīng)網(wǎng)絡(luò)是遞歸神經(jīng)網(wǎng)絡(luò)(RNN)的一種變體，LSTM擴展了其記憶能力[9]。圖1給出了LSTM內(nèi)部原理結(jié)構(gòu)示意圖，⊕代表兩個向量的加法運算，?代表兩個向量的點乘運算，σ為sigmoid激活函數(shù)，tanh為雙曲正切激活函數(shù)。

圖1 LSTM內(nèi)部原理

LSTM神經(jīng)網(wǎng)絡(luò)設(shè)計的核心是圍繞“記憶能力”開展的，為了實現(xiàn)這種記憶功能，LSTM引入了3個“門”的概念，即遺忘門ft、輸入門it和輸出門ot。遺忘門的主要功能是決定應(yīng)該保留或者拋棄哪些信息，輸入門用來更新單元狀態(tài)，輸出門可以控制下個狀態(tài)的隱藏信息的值。利用LSTM能有效捕捉長期時變信息的優(yōu)勢，以單顆衛(wèi)星前期鐘差時間序列為輸入，就能預(yù)測未來鐘差數(shù)據(jù)。

1.4 精度評價指標

本文以采用均方根誤差(RMS)、極差(Range)、小誤差概率(P參數(shù))、后驗方差比(C參數(shù))等4項指標來評價衛(wèi)星鐘差預(yù)報模型效果[1]；RMS能夠反映預(yù)報值與真值之間差異的分布狀態(tài)，Range可以反映出最大與最小預(yù)報誤差之差的絕對值大??；P參數(shù)與C參數(shù)可以通過對預(yù)報殘差考察進而精準判斷與公正評價預(yù)報模型的精度，表1給出了文獻[10]具體的精度等級評價標準，P參數(shù)越大則說明誤差較小的情況越多，C參數(shù)越小表示模型精度越高。

表1 精度等級評價標準

2 衛(wèi)星鐘差預(yù)報實驗分析

2.1 數(shù)據(jù)來源

本文選取武漢大學(xué)IGS數(shù)據(jù)中心(全球第5個國際GNSS服務(wù)組織授權(quán)負責建設(shè)的IGS數(shù)據(jù)中心)提供的2020年5月30日至2020年6月13日共15 d 33顆BDS衛(wèi)星精密鐘差產(chǎn)品(ftp://igs.ensg.ign.fr/pub/igs/products/mgex/)，表2給出了該33顆衛(wèi)星的星載原子鐘類型及衛(wèi)星編號，按照15 min采樣間隔對原始30 s采樣間隔鐘差序列進行提取，共獲取47 520個鐘差數(shù)據(jù)；短期鐘差預(yù)報思路是先讀取精密鐘差數(shù)據(jù)，利用文獻[11]衛(wèi)星鐘單天解基準處理方法消除天與天產(chǎn)品之間存在的系統(tǒng)性鐘差基準偏差，同時利用文獻[12]中改進的巴爾達法(Baarda)進行粗差和頻跳等異常值的探測、剔除與補充，然后建立PLS、DPLS、LSTM預(yù)報模型，對于每種預(yù)報模型下每顆衛(wèi)星，根據(jù)建模數(shù)據(jù)的不同分為兩類實驗：①均選擇第15天共24、96個數(shù)據(jù)為預(yù)報驗證數(shù)據(jù)，前14天距第15天最近的連續(xù)nd數(shù)據(jù)為建模數(shù)據(jù)，n從1～14變化，相應(yīng)預(yù)報實驗分別簡記為第1，2，…，14組實驗；②連續(xù)14 d 的短期預(yù)報(當前1天的數(shù)據(jù)為訓(xùn)練樣本，后1天的6 h、24 h數(shù)據(jù)為測試樣本)；最后按照式(7)、(8)、(11)、(12)進行精度評定。

表2 選取的BDS-2和BDS-3星載原子鐘類型及衛(wèi)星編號一覽表

2.2 實驗分析

2.2.1 建模數(shù)據(jù)量與預(yù)報精度之間關(guān)系分析

圖2為3種模型第1，2，…，14組實驗6 h，24 h預(yù)報結(jié)果中Range，RMS，P參數(shù)，C參數(shù)的變化情況。

從圖2可以看出：

圖2 PLS(藍)、DPLS(紅)、LSTM(綠)預(yù)報模型6 h、24 h預(yù)報結(jié)果

1)預(yù)報模型不同，建模數(shù)據(jù)量與預(yù)報精度之間的關(guān)系不盡相同；14組實驗中，6 h與24 h預(yù)報結(jié)果的變化規(guī)律基本相同；

2)PLS模型中，隨著建模數(shù)據(jù)量增加，預(yù)報精度變化不大，第8組精度相對最好，其他組精度變化呈現(xiàn)以第8組為對稱軸，逐漸向兩邊降低的趨勢，第一組精度相對較差，24 h預(yù)報結(jié)果中Range在4.63×10-3～7.06×10-3ns變化； RMS在2.56×10-3～3.10×10-3ns變化；P參數(shù)在0.98～1.00變化，C參數(shù)在2.88×10-6～4.38×10-6之間變化；24 h預(yù)報結(jié)果中每組實驗的精度等級均為一級，6 h預(yù)報結(jié)果中除1組實驗外其他組實驗的精度均為二級；

3)DPLS模型中，隨著建模數(shù)據(jù)量增加，預(yù)報精度變化較小，第1，3，8組精度相對較好，第2、4組精度相對較差，24 h預(yù)報結(jié)果中Range在9.85×10-4～1.11×10-4ns變化； RMS在1.90×10-3～2.21×10-3ns變化；P參數(shù)在0.99～1.00之間變化，C參數(shù)在1.92×10-6～2.45×10-6之間變化；24 h預(yù)報結(jié)果中每組實驗的精度等級均為一級，6 h預(yù)報結(jié)果中每組實驗的精度接近一級或為一級；

4)LSTM模型中，隨著建模數(shù)據(jù)量增加，預(yù)報精度呈現(xiàn)從低到高的變化趨勢，第11組預(yù)報精度相對最好，前6組的精度相對較差，第7～9組精度變化不大，第10～14組精度變化不大；24 h預(yù)報結(jié)果中Range均值可達到3.15×10-3～12.00×10-3ns；RMS均值達到1.70×10-3～6.37×10-3ns；P參數(shù)均值第10～14組均為1，其他組均近似為1，C參數(shù)均值可達到1.60×10-6～24.51×10-6，每組實驗的精度等級均為一級；6 h預(yù)報結(jié)果中除1～3組外，其他組實驗的精度接近為一級；

5)該3種模型均適合短期預(yù)報，PLS、DPLS模型并不需要較長的歷史數(shù)據(jù)就能獲得較好且穩(wěn)定的預(yù)報結(jié)果，主要原因是該算法能夠?qū)?shù)據(jù)進行主成分、相關(guān)性、線性回歸分析，能利用較少信息來歸納或解釋數(shù)據(jù)；DPLS模型較PLS精度更高，主要原因是通過對一次差序列進行改進的Baarda法粗差探測、剔除與補齊，不僅削弱了鐘差相位數(shù)據(jù)中趨勢項和粗差數(shù)據(jù)對模型訓(xùn)練的影響，還能減少輸入數(shù)據(jù)有效位數(shù)長度，使輸入的數(shù)據(jù)更能表征數(shù)據(jù)的變化特點，減少模型的冗余；LSTM模型需要較長的歷史數(shù)據(jù)才能達到最優(yōu)，超過相應(yīng)時長，LSTM模型精度變化微弱，主要原因是LSTM模型具有長短期記憶功能，能有效利用時間序列數(shù)據(jù)中長距離依賴信息的進行預(yù)測，但是在記憶及處理超過1 000量級的數(shù)據(jù)序列時會比較棘手，提升預(yù)報精度的能力受到限制。

2.2.2 預(yù)報耗時分析

圖3為3種模型第1，2，…，14組實驗24 h預(yù)報耗時變化情況，可以看出：

圖3 PLS(藍)、DPLS(紅)、LSTM(綠)預(yù)報模型24 h預(yù)報耗時情況

第14組實驗，PLS、DPLS、LSTM模型短期預(yù)報的耗時分別為0.24 s、0.84 s、198.65 min；LSTM模型耗時較長，隨著建模數(shù)據(jù)量增加，耗時逐步增長，最小、最大耗時分別為8.34 min、198.65 min，第3組實驗的耗時是第1組實驗的4.08倍，第14組是第1組的23.82倍數(shù)，主要原因是該模型無法實現(xiàn)并行優(yōu)化，訓(xùn)練速度較慢，時間跨度較大時，計算量會很大、很耗時；PLS、DPLS模型的耗時相當，最大耗時分別為0.25 s、0.90 s，最小耗時為0.21 s、0.76 s，均可以滿足s級預(yù)報要求，主要原因為該兩種模型可實現(xiàn)高維變量降維處理，具有較優(yōu)的實時性、魯棒性和泛化能力；可以看出，在滿足預(yù)報精度的情況下，選擇合適的預(yù)報方法、建模數(shù)據(jù)量是很有必要的。

2.2.3 衛(wèi)星鐘類型與預(yù)報精度之間關(guān)系分析

表3、表4分別為按照衛(wèi)星鐘類型，連續(xù)14 d的短期預(yù)報(當前1天數(shù)據(jù)為訓(xùn)練樣本，后1天6 h、24 h數(shù)據(jù)為測試樣本)中Range，RMS，P參數(shù)，C參數(shù)的統(tǒng)計結(jié)果；圖4為第1組實驗(以第14天數(shù)據(jù)為建模數(shù)據(jù)，第15天24 h數(shù)據(jù)為預(yù)報驗證數(shù)據(jù))得到的預(yù)報殘差序列圖。

表3 連續(xù)14 d短期預(yù)報Range，RMS統(tǒng)計結(jié)果

表4 連續(xù)14 d短期預(yù)報P參數(shù)，C參數(shù)統(tǒng)計結(jié)果

圖4 第一組實驗預(yù)報殘差序列圖

從表3、表4與圖4可以看出：

1)衛(wèi)星鐘類型不同，預(yù)報精度不同。6 h預(yù)報時長下，PLS、DPLS模型中BDS-3 Rb-II、BDS-3 PHM衛(wèi)星鐘預(yù)報精度相當，24 h預(yù)報時長下，BDS-3 PHM衛(wèi)星鐘預(yù)報精度略優(yōu)于BDS-3 Rb-II；兩種預(yù)報時長下，LSTM模型中BDS-3 PHM衛(wèi)星鐘預(yù)報精度明顯優(yōu)于BDS-3 Rb-II，相應(yīng)的RMS在6 h、24 h預(yù)報中分別提升了45.47%、56.09%；PLS、DPLS、LSTM模型中BDS-3衛(wèi)星預(yù)報精度明顯優(yōu)于BDS-2衛(wèi)星，6 h預(yù)報中RMS分別提升了63.47%、58.89%、68.10%，24 h預(yù)報中RMS分別提升了66.37%、66.91%、58.44%；每種模型預(yù)報結(jié)果中，BDS-2 Rb鐘較其他兩種類型的衛(wèi)星鐘，Range明顯較大，P參數(shù)明顯較小，C參數(shù)略大，說明其預(yù)報結(jié)果的穩(wěn)定性相對較差；

2)預(yù)報方案不同，預(yù)報精度不同。整體上DPLS預(yù)報精度最好，PLS次之，LSTM預(yù)報精度相對較差，24 h預(yù)報時長下，DPLS較PLS、LSTM的RMS分別提升了33.48%、42.68%，Range分別改善了29.77%、42.95%，P參數(shù)分別改善了12.76%、16.16%，C參數(shù)分別改善了42.20%、67.37%，6 h預(yù)報時長下，預(yù)報結(jié)果表現(xiàn)出了類似的精度改善情況；6 h預(yù)報中，DPLS較PLS，對于不同衛(wèi)星類型，Range與RMS的改善情況由高到低的順序為BDS-2 Rb、BDS-3 Rb-II、BDS-3 PHM，而在24 h預(yù)報中，不同衛(wèi)星類型的Range與RMS提升情況變化較?。?種預(yù)報時長下，DPLS較LSTM，6 h預(yù)報時長下，BDS-3 PHM鐘較其他兩種類型衛(wèi)星鐘，Range與RMS的改善程度明顯較小，在24 h預(yù)報時長下，BDS-3 Rb-II鐘較其他兩種類型衛(wèi)星鐘，Range與RMS的改善程度明顯較高。

從以上分析可以看出，DPLS模型的預(yù)報精度和穩(wěn)定性均優(yōu)于其他兩種模型，可以較好的應(yīng)用于超快速鐘差預(yù)報中，對實時精密單點定位具有重要意義。

3 結(jié) 論

本文基于武漢大學(xué)IGS數(shù)據(jù)中心連續(xù)15 d 33顆 事后BDS-2/BDS-3精密衛(wèi)星鐘差產(chǎn)品，利用DPLS算法進一步拓展了基于機器學(xué)習算法的衛(wèi)星鐘差短期預(yù)報，以PLS、LSTM模型為參照，分析了建模數(shù)據(jù)量、衛(wèi)星鐘類型與預(yù)報精度間的關(guān)系，建模數(shù)據(jù)量與預(yù)報耗時之間的關(guān)系，DPLS模型短期預(yù)報性能，得出以下結(jié)論：

1)該3種模型均適合短期預(yù)報；PLS、DPLS模型中，隨著建模數(shù)據(jù)量增加，預(yù)報精度變化不大，且不需要較長的歷史數(shù)據(jù)就能獲得較好且穩(wěn)定的預(yù)報結(jié)果，而LSTM模型預(yù)報精度呈現(xiàn)從低到高的變化趨勢，需要較長的歷史數(shù)據(jù)才能達到最優(yōu)，超過相應(yīng)時長，LSTM模型預(yù)報精度提升微弱；

2)PLS、DPLS模型耗時較少，LSTM模型耗時較長，利用連續(xù)14 d歷史數(shù)據(jù)進行短期預(yù)報的耗時分別為0.24 s、0.84 s、198.65 min；

3)利用一天的歷史數(shù)據(jù)進行短期預(yù)報的實驗表明，24 h預(yù)報時長下，DPLS較PLS、LSTM的RMS分別提升了33.48%、42.68%，Range分別改善了29.77%、42.95%，P參數(shù)分別改善了12.76%、16.16%，C參數(shù)分別改善了42.20%、67.37%，6h預(yù)報時長下，預(yù)報結(jié)果表現(xiàn)出了同一量級的精度改善情況；PLS、DPLS、LSTM模型中BDS-3衛(wèi)星預(yù)報精度明顯優(yōu)于BDS-2衛(wèi)星，6 h預(yù)報中RMS分別提升了63.47%、58.89%、68.10%，24 h預(yù)報中RMS分別提升了66.37%、66.91%、58.44%。

本研究結(jié)果對超快速鐘差產(chǎn)品高精度預(yù)報有一定參考意義；后期可基于PLS的改進算法(穩(wěn)健偏最小二乘、核偏小二乘、二階偏最小二乘)，從樣本的預(yù)處理(對樣本數(shù)據(jù)進行協(xié)方差的穩(wěn)健估計或奇異點的識別和剔除等處理)、特征選擇(過濾器模式，封裝器模式，嵌入式模式)等方面入手，來克服PLS或DPLS在鐘差預(yù)報中存在的過擬合問題和提高模型泛化能力，從而提升預(yù)報模型的精度；同時，可采用LSTM與Kalman模型相結(jié)合的鐘差預(yù)報方法來提高預(yù)報精度和穩(wěn)定性，即利用 LSTM架構(gòu)準確捕捉輸入輸出之間的狀態(tài)轉(zhuǎn)移關(guān)系，結(jié)合Kalman濾波的狀態(tài)空間模型對衛(wèi)星鐘差進行預(yù)測，既可以在一定程度上克服訓(xùn)練數(shù)據(jù)序列長度對傳統(tǒng) LSTM神經(jīng)網(wǎng)絡(luò)的影響，又能擺脫Kalman濾波對動力學(xué)模型的依賴使得模型可以從輸入數(shù)據(jù)中學(xué)習。