附有周期項的二次多項式LASSO鐘差預(yù)報模型

李坤，王潛心，閔揚(yáng)海，龔佑興，苗偉，程彤

(1.中國礦業(yè)大學(xué) 環(huán)境與測繪學(xué)院，江蘇 徐州 221116；2.蘇州市房地產(chǎn)市場與交易管理中心，江蘇 蘇州 215002；3.國防科學(xué)技術(shù)大學(xué) 指揮軍官基礎(chǔ)教育學(xué)院，湖南 長沙 410072)

0 引 言

北斗衛(wèi)星導(dǎo)航定位系統(tǒng)星載原子鐘因自身時頻特性的復(fù)雜性和外部因素影響，導(dǎo)致其系統(tǒng)時間與地面控制系統(tǒng)的原子鐘鐘面值存在偏差，該偏差即為衛(wèi)星鐘差[1-2]。目前，國際GNSS監(jiān)測評估系統(tǒng)(international GNSS monitoring &assessment system,IGMAS)最終解算的鐘差產(chǎn)品雖精度較高[3]，但解算具有一定滯后性,無法滿足實時或近實時用戶，因此，開發(fā)預(yù)報高精度鐘差的產(chǎn)品非常重要。

目前，關(guān)于鐘差預(yù)報的模型主要有二次多項式模型(quadratic polynomial,QP)、灰色模型(gray model,GM)、譜分析模型(spectrum analysis,SA)、時間序列模型(autoregressive integrated moving average,ARIMA)、卡爾曼濾波模型(kalman filtering,KF)[1,4]、神經(jīng)網(wǎng)絡(luò)模型(back propagation network,BP)等?；谛l(wèi)星鐘的物理特性，國內(nèi)外學(xué)者做了大量關(guān)于二次多項式預(yù)報模型的研究，并取得了豐碩成果[5-6]。鄭作亞等[7]、楊定江等[8]在二次多項式模型基礎(chǔ)上，增加了周期項因素，構(gòu)造附有周期項的二次多項式預(yù)報模型；王甫紅等[9]在附有周期項的二次多項式模型基礎(chǔ)上，提出了一種基于鐘差變化率擬合建模的衛(wèi)星鐘差預(yù)報方法；MAO Y等[10]在二次多項式模型基礎(chǔ)上，采用一種顧及衛(wèi)星間相關(guān)性的鐘差預(yù)報模型，并對模型中涉及的主要算法進(jìn)行了推導(dǎo)；HANG G等[11]、艾青松等[12]針對各天體間鐘差基準(zhǔn)偏差，在二次多項式模型基礎(chǔ)上提出了起點偏差修正法；于燁等[13]、王宇譜等[14]考慮了二次多項式預(yù)報模型隨機(jī)性變化的特點，對預(yù)報模型進(jìn)行了進(jìn)一步優(yōu)化。常規(guī)的二次多項式模型在求取模型系數(shù)時，常采用最小二乘估計(least squares estimation,LSQ)算法解算模型系數(shù)，但往往導(dǎo)致預(yù)報模型過擬合，因此，本文采用LASSO算法進(jìn)行模型系數(shù)的解算。先分別使用LSQ和LASSO算法求得模型參數(shù)，比較兩種模型擬合殘差，分析模型擬合精度，求取模型周期項，然后對附有周期項的各模型進(jìn)行預(yù)報，比較6，12，18，24 h預(yù)報精度。

1 二次多項式模型參數(shù)求解方法

1.1 LSQ算法

不考慮殘差時，可表示為

y=a0+a1t+a2t2，

(1)

式中：a0為鐘差；a1為鐘速或頻偏；a2為頻漂；t為時間。

二次多項式預(yù)報模型的誤差方程可表示為

v=Bx-l，

(2)

式中：v為誤差項；B為鐘差歷元系數(shù)矩陣；l為已知鐘差；x為待求模型參數(shù)。

針對上述預(yù)報模型和誤差方程，根據(jù)最小二乘原理求得模型系數(shù)為[15]

x=(BTPB)-1BTPl，

(3)

其中，P為權(quán)陣(本文所取得的權(quán)陣為單位陣)。

1.2 LASSO算法

LASSO算法是T.Robert[16]提出的一種處理復(fù)共線性數(shù)據(jù)的有偏估計。LASSO算法的代價函數(shù)可表示為

J(x)=(Y-Ax)2+λ‖x‖1,

(4)

式中：A為鐘差歷元系數(shù)矩陣；Y為已知鐘差；λ為懲罰系數(shù)；‖x‖1為x的1范數(shù)。

對式(4)右邊第一項求偏導(dǎo),得

-2mj+2xjnj，

(5)

因式(4)第二項不可導(dǎo)，對其使用次導(dǎo)數(shù)，

令式(5)，式(6)相加等于0，則

(7)

由式(7)可得

(8)

綜合考慮懲罰系數(shù)的有效性和算法的運(yùn)行速率，選取經(jīng)驗系數(shù)0.01,0.1,1,2,3,5,10,100,200,500為懲罰系數(shù)。為了達(dá)到理論最優(yōu)，懲罰系數(shù)應(yīng)越大，相應(yīng)附加懲罰系數(shù)的參數(shù)越小，甚至為0，從而達(dá)到減少變量和降低方程維度的目的。相關(guān)學(xué)者研究了窮舉法選取最優(yōu)懲罰系數(shù)，但耗時較長，本文選取常用經(jīng)驗系數(shù)，通過模型迭代和篩選，獲得最佳懲罰系數(shù)，并將懲罰系數(shù)回代到模型中進(jìn)行計算。

2 二次多項式鐘差預(yù)報策略

2.1 數(shù)據(jù)預(yù)處理

GNSS原始鐘差數(shù)據(jù)存在數(shù)據(jù)異常問題，同時由于外部因素影響，如設(shè)備老化、外部環(huán)境變化等導(dǎo)致原始鐘差數(shù)據(jù)中存在粗差和鐘跳等情況，對后續(xù)鐘差預(yù)報和分析帶來影響[18]。針對上述情況，需要對鐘差進(jìn)行預(yù)處理，得到有效數(shù)據(jù)，為后續(xù)數(shù)據(jù)的正確預(yù)報和分析提供可靠基礎(chǔ)。

鐘差數(shù)據(jù)預(yù)處理方法包括中位數(shù)法(median absolute deviation,MAD)和Baarda粗差探測法。中位數(shù)法簡單、方便，但對粗差的大小不敏感，不能探測出小粗差。Baarda粗差探測法通過對驗后殘差進(jìn)行處理，能有效探測出小粗差，但在進(jìn)行鐘差系數(shù)擬合求取鐘差殘差時，常因鐘跳等異常使得鐘差值無法準(zhǔn)確獲得擬合系數(shù)[19]。

針對Baarda粗差探測法和中位數(shù)法的優(yōu)缺點，對Baarda算法進(jìn)行優(yōu)化和改進(jìn)，先將鐘差相位數(shù)據(jù)轉(zhuǎn)化為頻率數(shù)據(jù)，具體數(shù)學(xué)模型為

(9)

式中：x(t)，y(t)分別為相位和頻率數(shù)據(jù)；εt和ωt分別為相位和頻率數(shù)據(jù)的時間序列擬合殘差。

對頻率數(shù)據(jù)進(jìn)行最小二乘平差估計得到模型的擬合殘差，計算單位權(quán)方差，判斷最大殘差是否超過閾值。實際運(yùn)算中有部分殘差超過閾值，但最小二乘算法具有殘差平均功能，不能剔除超限的其他殘差，而是重新組建方程，然后再次進(jìn)行殘差的剔除。

最后對剔除殘差的鐘差頻率數(shù)據(jù)進(jìn)行二次多項式擬合，補(bǔ)全缺失的數(shù)據(jù)，同時將頻率數(shù)據(jù)恢復(fù)為相位數(shù)據(jù)，便于后續(xù)鐘差預(yù)報等工作。圖1為鐘差數(shù)據(jù)預(yù)處理流程。

圖1 鐘差數(shù)據(jù)預(yù)處理流程Fig.1 clock offset data preprocessing

為檢驗改進(jìn)Baarda算法的有效性和可行性，本文采用德國地學(xué)中心(GFZ)發(fā)布的GBM鐘差數(shù)據(jù)，選取2017年10月1日，即gbm19690.clk鐘差數(shù)據(jù)，采樣間隔為30 s，進(jìn)行鐘差數(shù)據(jù)預(yù)處理。如圖2所示，為了準(zhǔn)確分析和直觀顯示改進(jìn)的Baarda算法對數(shù)據(jù)處理的效果，對GBM鐘差數(shù)據(jù)中的C02號衛(wèi)星人工添加粗差數(shù)據(jù)，利用Barrda算法進(jìn)行粗差探測，并修復(fù)鐘差數(shù)據(jù)中的粗差和鐘跳數(shù)，對于頻漂大的鐘差數(shù)據(jù)，很難通過常規(guī)方法進(jìn)行預(yù)處理，而通過圖2可以看出，改進(jìn)后的Baarda算法能夠很好地識別和處理粗差剔除數(shù)據(jù)。

圖2 C02衛(wèi)星頻率數(shù)據(jù)預(yù)處理效果Fig.2 Frequency data preprocessing result of C02 satellite

2.2 擬合殘差分析

自鄭作亞等[7]分析了附有周期項的二次多項式預(yù)報模型以來，周期項誤差已成為鐘差基本屬性之一，鐘差主要分為趨勢項和周期項兩部分，鑒于此，本文采用附有周期項的二次多項式預(yù)報模型，具體模型為

(10)

式中：P為周期項個數(shù)；Ai，Bi，t′分別為正、余弦函數(shù)的振幅和周期項；εi為附加周期項后的擬合殘差。

采用IGMAS 2018年6月18日快速鐘差產(chǎn)品進(jìn)行模型參數(shù)求解。首先對鐘差數(shù)據(jù)預(yù)處理，獲得有效的鐘差數(shù)據(jù)，然后分別采用LSQ和LASSO算法求得附有周期項的二次多項式預(yù)報模型系數(shù)，用求解的模型參數(shù)分別擬合當(dāng)天鐘差，并與當(dāng)天實際鐘差作差，求解模型擬合精度。圖3～4為各取兩顆GEO衛(wèi)星(C01，C04)，IGSO衛(wèi)星(C07，C09)，MEO衛(wèi)星(C12,C14)進(jìn)行模型精度分析。

圖3 模型擬合殘差時間序列Fig.3 Fitting residual time series diagram

由圖3可以看出，LSQ算法求得的自擬合殘差分布在[-2,2]，且相對集中，而LASSO算法求得的自擬合殘差則明顯分散，且大于LSQ算法的?？梢缘贸?，LSQ算法求取的模型系數(shù)自擬合程度明顯優(yōu)于LASSO算法的。

為了定性分析LSQ和LASSO算法的自擬合精度，分別統(tǒng)計圖3衛(wèi)星精度情況，并繪制模型擬合精度直方圖(圖4)。由圖4可以得出，LSQ算法擬合精度整體優(yōu)于LASSO算法的，除C07衛(wèi)星外，LSQ算法求解的模型擬合精度較高，小于0.5 ns，可以看出，使用LSQ算法進(jìn)行模型擬合，求取的模型具有很高的擬合精度。但是否存在過擬合導(dǎo)致模型錯誤，從而使得模型在預(yù)報時偏離太大，下文將對該現(xiàn)象進(jìn)行具體分析。

圖4 模型擬合精度Fig.4 Model fitting accuracy

3 實驗結(jié)果分析

為驗證LASSO算法求得預(yù)報模型的可行性，采用IGMAS 2018年6月18日的快速鐘差產(chǎn)品進(jìn)行模型系數(shù)求解，以2018年6月19日的快速鐘差產(chǎn)品為真值，進(jìn)行預(yù)報精度的分析。

首先對2018年6月18日的鐘差數(shù)據(jù)進(jìn)行預(yù)處理，獲取有效頻率數(shù)據(jù)，然后對鐘差數(shù)據(jù)進(jìn)行附有周期項的二次多項式建模，用LSQ和LASSO算法分別求解模型參數(shù)，用獲得的模型求解鐘差預(yù)報產(chǎn)品，對求取的鐘差預(yù)報產(chǎn)品與6月19日的鐘差產(chǎn)品進(jìn)行對比分析，獲取鐘差預(yù)報精度數(shù)據(jù)。最后聯(lián)合擬合殘差數(shù)據(jù)和鐘差預(yù)報產(chǎn)品精度數(shù)據(jù)進(jìn)行過擬合分析，分析LSQ算法是否存在過擬合導(dǎo)致模型錯誤。實驗分析流程見圖5。

圖5 實驗分析流程圖Fig.5 Flow chart of experimental analysis

統(tǒng)計分析上述附有周期項的LSQ和LASSO算法的鐘差預(yù)報模型中預(yù)報的鐘差，用于后續(xù)過擬合現(xiàn)象的具體分析。

根據(jù)實際情況，分別統(tǒng)計各衛(wèi)星6，12，18，24 h鐘差預(yù)報精度。圖6為兩種方法求解的模型預(yù)報殘差時間序列，為更直觀地分析兩種算法的預(yù)報結(jié)果，分別統(tǒng)計兩種方法的預(yù)報均方根誤差(root mean square，RMS)，并繪制直方圖(圖7)。兩種算法不同時間段的鐘差預(yù)報精度數(shù)據(jù)如表1所示。

表1 兩種方法預(yù)報精度統(tǒng)計Tab.1 Prediction accuracy statistics of the two methods

由圖6可知，LSQ和LASSO算法求解的鐘差預(yù)報殘差為[-40，40]ns，LASSO算法求解的模型預(yù)報發(fā)散程度明顯小于LSQ算法的。由圖7可知，LASSO算法求解的預(yù)報模型，各衛(wèi)星的預(yù)報精度都優(yōu)于LSQ算法的。根據(jù)圖3～4和圖6～7對比分析可得，LSQ算法求解的模型擬合精度高，但預(yù)報精度差。

圖6 鐘差預(yù)報殘差時間序列圖Fig.6 Time series diagram of residuals of clock offset prediction

圖7 鐘差預(yù)報精度均方根誤差直方圖Fig.7 RMS histograms of clock offset prediction accuracy

由于本文采用附有周期項的二次多項式預(yù)報模型，周期項數(shù)據(jù)未做篩選，導(dǎo)致模型變量非常多。常規(guī)無偏LSQ算法在處理多變量數(shù)據(jù)時易產(chǎn)生過擬合現(xiàn)象，即模型本身擬合精度非常高，但將模型外推時其精度很差，因而不能很好地表現(xiàn)模型的一般性。LASSO算法的最大優(yōu)勢為添加了懲罰項，可以降低模型復(fù)雜度，即并非所有變量都放入模型中，而是將部分參數(shù)的系數(shù)變?yōu)?，有選擇地將一些性能更好的參數(shù)加入模型中，從而降低模型維度，避免過擬合現(xiàn)象。

比較圖3、圖6可以看出，考慮所有周期項變量的LSQ算法自擬合精度非常高，是因為LSQ算法本身對所有周期進(jìn)行建模，相比LASSO算法剔除一些變量而言，LSQ算法自擬合精度顯著高于LASSO算法的，但數(shù)據(jù)變量龐大和觀測數(shù)據(jù)有限之間的矛盾導(dǎo)致LSQ算法求解預(yù)報模型時存在過擬合現(xiàn)象，造成求解的模型預(yù)報精度不高，而LASSO算法能很好地避免這一現(xiàn)象發(fā)生。

由表1可得，LASSO算法求解的模型預(yù)報精度提高較顯著，除預(yù)報時長6 h的C09，C12衛(wèi)星和預(yù)報時長12 h的C12衛(wèi)星外，其余的提升效率均優(yōu)于10%，且隨著預(yù)報時間增加，LASSO算法求解的模型預(yù)報精度提升效果越顯著。

4 結(jié) 論

(1)LSQ算法求解的模型系數(shù)自擬合精度顯著優(yōu)于LASSO算法的，但模型預(yù)報精度在預(yù)報時長6，12，18，24 h時均低于LASSO算法的，因此，使用LSQ算法進(jìn)行模型參數(shù)求解，存在過擬合現(xiàn)象，從而導(dǎo)致預(yù)報精度不高。

(2)LASSO算法求解的模型預(yù)報精度較LSQ算法提升明顯，除預(yù)報時長6 h的C09，C12衛(wèi)星和預(yù)報時長12 h的C12衛(wèi)星外，其余均提升超過10%。從總體趨勢上看，隨著預(yù)報時間增加，LASSO算法求解的模型優(yōu)勢更為明顯，其預(yù)報發(fā)散程度顯著小于LSQ算法解算的模型。因此，LASSO算法可以抑制過擬合現(xiàn)象，大幅提升預(yù)報模型的準(zhǔn)確性，對鐘差預(yù)報的改進(jìn)具有顯著意義。