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

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

0 引 言

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

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

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

1.1 LSQ算法

不考慮殘差時(shí)，可表示為

y=a0+a1t+a2t2，

(1)

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

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

v=Bx-l，

(2)

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

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

x=(BTPB)-1BTPl，

(3)

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

1.2 LASSO算法

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

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

(4)

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

對(duì)式(4)右邊第一項(xiàng)求偏導(dǎo),得

-2mj+2xjnj，

(5)

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

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

(7)

由式(7)可得

(8)

綜合考慮懲罰系數(shù)的有效性和算法的運(yùn)行速率，選取經(jīng)驗(yàn)系數(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ù)，但耗時(shí)較長(zhǎng)，本文選取常用經(jīng)驗(yàn)系數(shù)，通過模型迭代和篩選，獲得最佳懲罰系數(shù)，并將懲罰系數(shù)回代到模型中進(jìn)行計(jì)算。

2 二次多項(xiàng)式鐘差預(yù)報(bào)策略

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

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

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

針對(duì)Baarda粗差探測(cè)法和中位數(shù)法的優(yōu)缺點(diǎn)，對(duì)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í)間序列擬合殘差。

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

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

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

為檢驗(yàn)改進(jìn)Baarda算法的有效性和可行性，本文采用德國(guó)地學(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算法對(duì)數(shù)據(jù)處理的效果，對(duì)GBM鐘差數(shù)據(jù)中的C02號(hào)衛(wèi)星人工添加粗差數(shù)據(jù)，利用Barrda算法進(jìn)行粗差探測(cè)，并修復(fù)鐘差數(shù)據(jù)中的粗差和鐘跳數(shù)，對(duì)于頻漂大的鐘差數(shù)據(jù)，很難通過常規(guī)方法進(jìn)行預(yù)處理，而通過圖2可以看出，改進(jìn)后的Baarda算法能夠很好地識(shí)別和處理粗差剔除數(shù)據(jù)。

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

2.2 擬合殘差分析

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

(10)

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

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

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

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

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

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

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

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

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

圖5 實(shí)驗(yàn)分析流程圖Fig.5 Flow chart of experimental analysis

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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