基于BDS-2/BDS-3聯(lián)合處理的北斗超快速鐘差預報優(yōu)化策略

馬昌忠，王潛心，胡 超，閔揚海，王澤杰

(1. 中國礦業(yè)大學自然資源部國土環(huán)境與災害監(jiān)測重點實驗室，徐州 221116;2. 中國礦業(yè)大學環(huán)境與測繪學院，徐州 221116)

0 引言

目前，國際GNSS監(jiān)測與評估系統(tǒng)(interna-tional GNSS Monitoring and Assessment System, iGMAS)發(fā)布的超快速鐘差產(chǎn)品中，中國北斗衛(wèi)星導航系統(tǒng)(BeiDou Navigation Satellite System, BDS)觀測和預報部分精度分別約為0.60ns和6.00ns，遠遠不能滿足實時cm級應(yīng)用需求[1]。因此，在提高全球?qū)Ш叫l(wèi)星系統(tǒng)產(chǎn)品可用性、穩(wěn)定性和可靠性的前提下，有必要進一步完善BDS超快速衛(wèi)星鐘差產(chǎn)品質(zhì)量，尤其是預報鐘差精度。

當前，針對衛(wèi)星鐘差預報的研究主要集中在3個方面：1)衛(wèi)星鐘差序列預處理策略[2-4]；2)預報模型精化處理[5-7]；3)環(huán)境因素對鐘差建模的影響分析[8-9]。此外，在鐘差預報算法和策略研究中，學者主要關(guān)注鐘差序列長期預報模型研究，如擴展狀態(tài)模型[10]和人工神經(jīng)網(wǎng)絡(luò)預報[11]；多個短周期項疊加預報，如改進的迭代法以及基于恒星日濾波的單天內(nèi)鐘差變化量預報等[5]。通常衛(wèi)星鐘差預報模型采用趨勢項(多項式)與周期項表示[5,12]。同時，毛亞等提出了一種改進的預報策略[13]，即考慮星間相關(guān)性對模型參數(shù)估計的影響，一定程度上提升了BDS-2鐘差預報精度。

然而，BDS超快速鐘差建模中廣泛采用組合趨勢項與周期項的函數(shù)模型[2,14]，主要以全球定位系統(tǒng)(Global Positioning System，GPS)相關(guān)研究為參考。此外，黃觀文等提出了一種改進的BDS超快速鐘差預報策略，其中非線性項通過BP神經(jīng)網(wǎng)絡(luò)(Back Propagation Neural Network，BPNN)進行建模處理[1]。陳金平等評估了基于星間鏈路觀測數(shù)據(jù)的BDS-3試驗星鐘差預報特性[15]。此外，研究表明BDS-3配備的星載原子鐘頻率穩(wěn)定性較BDS-2提升了1個量級以上。胡超等通過引入星間相關(guān)性[16]，聯(lián)合預報了BDS-2和BDS-3超快速產(chǎn)品，并評估了星間相關(guān)性對超快速鐘差預報產(chǎn)品精度的影響[17]；基于1個月的實驗數(shù)據(jù)表明，引入星間相關(guān)性可實現(xiàn)BDS-2和BDS-3衛(wèi)星18h預報鐘差精度分別提高30.7%～47.3%和49.9%～59.3%。因此，為進一步提高北斗超快速預報鐘差產(chǎn)品質(zhì)量，后續(xù)研究中有必要對BDS-2和BDS-3鐘差聯(lián)合處理中的相關(guān)性系數(shù)進行深入挖掘分析與處理。

與當前廣泛使用的鐘差預報模型[1,13]相比，北斗衛(wèi)星鐘差建模中，模型參數(shù)估計存在2個主要問題需要優(yōu)化：1)建模過程中，通常將趨勢項與周期項分為兩部分進行獨立處理，不可避免地忽略了未知參數(shù)之間的相關(guān)性[16]；2)在建模中，只考慮顯著周期項(如BDS-2 2個周期、BDS-3 3個周期)[17]，將導致模型殘差無法建模。因此，為了獲得更準確且可靠的北斗鐘差預報模型，需要對BDS鐘差建模策略進行優(yōu)化處理，如精化鐘差模型項以及模型系數(shù)等。本文將利用一步求解策略，實現(xiàn)北斗各衛(wèi)星鐘差模型的統(tǒng)一估計；同時在一步模型參數(shù)估計中，對預報模型的精度、參數(shù)數(shù)目和參數(shù)類型(周期項和趨勢項個數(shù))進行權(quán)衡；為獲得最優(yōu)鐘差模型，模型參數(shù)估計中借助機器學習中稀疏建模的思想進行優(yōu)化處理[18]。

本文主要針對傳統(tǒng)鐘差預報模型建模策略的缺點，采用BDS-2/BDS-3鐘差序列聯(lián)合處理的方法，優(yōu)化北斗超快速衛(wèi)星鐘差產(chǎn)品。首先，基于稀疏建模方法提出了一種模型選擇方法，對傳統(tǒng)的鐘差預報兩步策略進行改進；其次，設(shè)計了一種BDS-2/BDS-3聯(lián)合處理策略，以充分利用BDS-2與BDS-3星間相關(guān)性實現(xiàn)參數(shù)解算增強；然后，考慮預報模型殘差序列的時空相關(guān)性，利用半變異函數(shù)構(gòu)建經(jīng)驗模型，精化北斗鐘差模型參數(shù)；最后，通過大量實驗對本文的優(yōu)化策略進行驗證。

1 北斗超快速鐘差預報優(yōu)化策略

1.1 BDS-2/BDS-3聯(lián)合處理模型

在對北斗原始鐘差序列進行嚴格預處理的基礎(chǔ)上，如粗差探測與修復和降噪等[16]，可將超快速解算鐘差序列設(shè)為L，則鐘差模型可表示為[1,13,16]

(1)

式中，ti和k分別表示第ti個歷元時間和第k顆衛(wèi)星；a0、a1和a2為趨勢項(多項式)系數(shù)；n為周期項總數(shù)，j為第j個周期；Aj、Bj和Tj表示周期項振幅及其相應(yīng)周期大小；ε(ti)表示模型殘差項。

式(1)中，周期項通常采用快速傅立葉變換算法計算得到。因此，根據(jù)得到的周期項，相應(yīng)的鐘差模型誤差方程可表示為

(2)

式中，Lk為鐘差減去周期項結(jié)果。需要注意的是，傳統(tǒng)鐘差序列建模策略中[16-17]，解算鐘差模型參數(shù)(趨勢項與周期項)分為2個獨立處理過程分別進行平差估計；理論上單個模型的所有參數(shù)需經(jīng)統(tǒng)一解算得到。同時，由于星載原子鐘的特性存在差異，所有鐘差模型參數(shù)的個數(shù)和類型不可能表示為統(tǒng)一形式；因此，為獲得更精確的鐘差模型，有必要細化北斗鐘差模型參數(shù)。

首先，一步估計式(2)中所有模型參數(shù)。為提高計算效率，同時解算所有衛(wèi)星多項式參數(shù)，假設(shè)有b顆衛(wèi)星、s個歷元，式(1)矩陣形式為

(3)

(4)

式(4)中，權(quán)矩陣P(t1)可表示為

(5)

rkb=

(6)

在鐘差模型參數(shù)估計中，有必要對權(quán)陣進行調(diào)整，以提高預報模型的精度，本文在權(quán)陣中引入時間變量，即

(7)

式中，Δt表示相鄰歷元的時間間隔。式(2)～式(7)，即為基于BDS-2/BDS-3星間相關(guān)性的北斗鐘差預報模型一步建模策略。

1.2 基于稀疏建模的鐘差模型估計

在1.1節(jié)的基礎(chǔ)上，鐘差模型參數(shù)估計考慮了衛(wèi)星間相關(guān)性，間接地優(yōu)化了模型參數(shù)求解過程。然而，參數(shù)估計中仍存在2個關(guān)鍵問題會導致待估參數(shù)無法獲得最優(yōu)解：1)由于建模中觀測值(鐘差序列)有限和模型過度擬合等因素影響，模型參數(shù)的穩(wěn)定性和精度不可避免地受到限制；2)不同類型星載原子鐘呈現(xiàn)出不同的特性(鐘差模型)[19]，無法基于統(tǒng)一的模型進行描述(相同的待估參數(shù))。為提高鐘差建模精度與穩(wěn)定性，有必要進一步優(yōu)化北斗衛(wèi)星鐘差建模策略。

本節(jié)利用機器學習中稀疏建模的思想進行鐘差模型精確構(gòu)建。理論上，引入所有的自變量進行建模處理可有效提高建模精度，但同時也增加了模型復雜性并導致過度擬合。在實際應(yīng)用中，為了提高模型的可靠性與預報精度，必須考慮模型自變量的合理選取，即模型選擇(或變量選擇)[20-21]。通常，模型選擇只能在有限的備選模型范圍內(nèi)進行，而稀疏建模可以有效地選擇建模數(shù)據(jù)與自變量(模型項)，從而實現(xiàn)模型預報精度和計算效率的提升。為對北斗鐘差序列進行稀疏建模處理，本文采用Tibshirani提出的最小絕對值收斂和選擇算子(Least absolute shrinkage and selection operator，Lasso)算法[22]。在Lasso算法中，通過最小化多變量線性回歸模型的損失函數(shù)，對回歸系數(shù)絕對值施加約束，實現(xiàn)了參數(shù)估計和模型選擇同步進行。根據(jù)Lasso算法原理，可實現(xiàn)對部分變量進行壓縮和消除；Lasso算法的具體實現(xiàn)是通過L1范數(shù)的最小二乘法，即

(8)

(9)

式中，λ為正則化系數(shù)。等式右邊相當于模型殘差平方和與正則化函數(shù)的組合；區(qū)別于Tikhonov正則化的L2范數(shù)，L1正則化范數(shù)通過自動模型選擇達到稀疏性的同時，不失函數(shù)模型的凸性(產(chǎn)生少量特征)。為求解式(9)，可通過迭代加權(quán)最小二乘法求解Lasso問題的數(shù)值解[23]，即式(9)可表示為

(10)

(11)

式中，U(m)是權(quán)陣(對角陣)，其元素根據(jù)第m-1次迭代解確定。因此

(12)

通過矩陣求逆公式，式(12)的解為

(13)

1.3 基于半變異函數(shù)的模型殘差處理

通過上述基于稀疏建模思想的一步求解鐘差模型參數(shù)，可較好地實現(xiàn)北斗鐘差預報模型構(gòu)建。但是受模型殘差以及隨機噪聲影響，建模過程中需對模型殘差序列進行深入分析。前期研究中，針對預報模型殘差序列提出了幾種改進的處理策略[1,13]；此外，學者組合了偏最小二乘法(Partial Least Square，PLS)和BPNN算法，并將其應(yīng)用于超快速鐘差預報中[14,17]；但缺少衛(wèi)星鐘差精度的相關(guān)信息(方差-協(xié)方差陣)，殘差序列建模及其預報難以達到預期精度[17]。本文從鐘差序列特點的角度考慮，即所有衛(wèi)星建模殘差序列不可避免地呈現(xiàn)時空相關(guān)性；因此，為充分挖掘BDS-2與BDS-3鐘差序列各衛(wèi)星潛在的時空特性，本節(jié)提出了利用半變異函數(shù)進行殘差序列優(yōu)化處理。有關(guān)半變異函數(shù)的更多細節(jié)已在相關(guān)研究中進行了詳細討論與總結(jié)[17,24-25]。

在實際建模處理中，通常使用改進的殘差半變異函數(shù)模型[19,24]

(14)

2γ(ti,tl)=2(C(0)-C(ti-tl))

(15)

式中，C表示參數(shù)協(xié)變差運算。C(tl-ti)可以通過式(16)計算得到

C(ti-tl)=Cov(ε(ti),ε(tl))

(16)

其中，Cov表示協(xié)方差運算；模型殘差的時間相關(guān)性將由半變異函數(shù)進行數(shù)值量化。

在得到每顆衛(wèi)星量化的殘差序列時間相關(guān)性之后，可擬合出經(jīng)驗半變異函數(shù)模型，如球形和指數(shù)型等?；诘玫降慕?jīng)驗半變異函數(shù)模型(如式(17)與式(18))，當h→∞時，2γ(h)→2C(0)；由式(15)可知，協(xié)變差C(h)可由經(jīng)驗函數(shù)模型確定。因此，在式(16)的基礎(chǔ)上可進一步構(gòu)造參數(shù)估計的權(quán)陣，并將其代入式(3)中，從而實現(xiàn)模型系數(shù)估計的精化處理。

綜上，本文從星間相關(guān)性、稀疏建模以及半變異函數(shù)等3個方面對北斗超快速預報鐘差模型進行了精化處理。首先，基于稀疏建模對鐘差模型參數(shù)進行一步估計，實現(xiàn)了模型的自動選擇，并改進了傳統(tǒng)鐘差分步建模策略；其次，考慮BDS-3星載原子鐘具有更強性能，通過提取星間相關(guān)性實現(xiàn)了模型參數(shù)的解算增強；最后，針對模型殘差序列，利用半變異函數(shù)對參數(shù)估計隨機模型進行了精化處理。

2 BDS-2/BDS-3衛(wèi)星鐘差預報實驗分析

本節(jié)將對改進的BDS-2/BDS-3衛(wèi)星超快速鐘差聯(lián)合預報策略進行實驗分析。實驗前，首先基于連續(xù)1個月(2019年第41～70天)的BDS-2和BDS-3觀測數(shù)據(jù)進行超快速聯(lián)合定軌實驗，得到了解算的超快速衛(wèi)星鐘差產(chǎn)品。鐘差預報實驗中，根據(jù)本文提出的改進策略，將實驗分為3組，分別從星間相關(guān)性、稀疏建模和半變異函數(shù)模型的角度對鐘差預報模型進行驗證分析?？紤]到北斗地球靜止軌道(Geostationary Earth Orbits，GEO)衛(wèi)星鐘差精度較低的特點，在實驗中剔除了相應(yīng)部分。

第1組實驗中包括4個對比方案，驗證了基于稀疏建模的鐘差序列一步建模策略：

方案1：基于趨勢項和周期項，建立鐘差預報模型，該方案中，趨勢項系數(shù)通過擬合鐘差序列得到，而周期項則通過殘差的快速傅里葉變換確定；

方案2：與方案1類似，其中BDS-2衛(wèi)星周期項參考黃觀文等[1]以及毛亞等[13]的研究結(jié)論，BDS-3則利用3個主要周期項建立鐘差預報模型；

方案3：與方案1和方案2相比，將快速傅里葉變換獲得的所有周期項加入鐘差序列建模中；

方案4：在方案3的基礎(chǔ)上，采用Lasso算法對模型參數(shù)進行一步估計，以驗證稀疏建模的有效性。

根據(jù)上述四種方案分別進行了1天弧長的BDS鐘差預報實驗，并以武漢大學iGMAS分析中心發(fā)布的BDS快速鐘差為參考進行精度分析。各方案中BDS-2與BDS-3衛(wèi)星的鐘差建模殘差以及相應(yīng)預報精度如圖1和圖2所示。

通過上述4套對比實驗方案可知：1)BDS-3衛(wèi)星鐘差建模殘差明顯小于BDS-2衛(wèi)星，間接說明了BDS-3星載原子鐘性能優(yōu)于BDS-2；2)方案1與方案2的擬合殘差近似相等，而方案3中利用所有周期項建模得到的模型殘差最??；3)基于圖2中鐘差預報精度的統(tǒng)計結(jié)果，對BDS-2而言，方案3較傳統(tǒng)方法(方案1與方案2)預報精度有所提高，且方案4(一步建模)略優(yōu)于方案3(除C13外)；對于BDS-3衛(wèi)星，方案3一步估計所得到的預報精度有所降低，但在稀疏建模之后(方案4)，在方案1與方案2的基礎(chǔ)上，一定程度實現(xiàn)了預報精度的提高。

圖1 年積日41天(2019)不同北斗衛(wèi)星鐘差擬合殘差Fig.1 Fitting residuals of BDS clock offsets for different satellites on DOY 41

圖2 四種方案連續(xù)30天北斗鐘差預報精度(18h)平均值Fig.2 Average daily RMS values of BDS predicted clock offsets for different satellites in 30-day experiments with four schemes (18 hours)

第2組實驗包括2個對比方案，驗證了引入星間相關(guān)性求解模型參數(shù)的可行性：

方案5：類似方案1，但在模型參數(shù)求解中加入了星間相關(guān)性系數(shù)，該方案具體實驗細節(jié)可參考作者前期研究成果[16-17]；

方案6：同樣地，在方案4的基礎(chǔ)上，將獲取的星間相關(guān)性系數(shù)加入衛(wèi)星鐘差模型參數(shù)估計中。

圖3所示為連續(xù)10天C14與C24衛(wèi)星預報鐘差的精度；表1則給出了1個月鐘差精度平均值(18h)及其精度提升率。實驗結(jié)果表明，星間相關(guān)性可顯著改善北斗衛(wèi)星預報鐘差精度；且與傳統(tǒng)建模方法(方案1)相比，BDS-2與BDS-3衛(wèi)星鐘差預報精度分別提升了32.9%和16.9%；而方案4考慮了星間相關(guān)性后，BDS-2與BDS-3衛(wèi)星鐘差18h預報精度分別提升了27.2%與28.6%。

圖3 不同方案連續(xù)10天鐘差預報精度Fig.3 Prediction accuracy of clock offsets in 10 consecutive days for different schemes

表1 不同方案北斗衛(wèi)星鐘差預報精度(ns)及其提升率(18h)

通過上述引入星間相關(guān)性與稀疏建模處理后，模型殘差中仍包含了鐘差序列的有效成分，其一定程度上降低了北斗鐘差預報的精度。第3組實驗主要分析了本文提出的超快速鐘差預報策略的正確性，并驗證了利用半變異函數(shù)精化模型參數(shù)估計的隨機模型的可行性。在鐘差預報實驗前，首先基于方案4中的模型殘差序列進行時空相關(guān)性估計，結(jié)果如圖4所示，其中虛線表示基于1個月殘差序列計算的實驗半變異函數(shù)值；其次，構(gòu)建了一個經(jīng)驗模型，即半變異函數(shù)模型，實驗中選擇球面模型進行殘差半變異函數(shù)模型構(gòu)建，即如式(17)與(18)所示

γ(h)BDS-2=

(17)

γ(h)BDS-3=

(18)

圖4 BDS-2和BDS-3方案4鐘差殘差序列半變異函數(shù)Fig.4 Variogram from the satellite clock offsets residuals of scheme 4 for BDS-2 and BDS-3 satellites

方案7：基于BDS-2與BDS-3聯(lián)合超快速鐘差序列，考慮星間相關(guān)性，通過方案1進行鐘差預報；

方案8：基于BDS-2與BDS-3聯(lián)合超快速鐘差序列，考慮星間相關(guān)性，通過方案2進行鐘差預報；

方案9：基于BDS-2與BDS-3聯(lián)合超快速鐘差產(chǎn)品，通過方案6進行鐘差預報，每個預報弧段長度為18h，并利用BPNN算法對模型殘差進行處理與預報，最后將兩部分預報鐘差序列合并；

方案10：類似于方案9，針對模型殘差部分采用灰色模型進行預報處理；

方案11：基于方案6，采用PLS+BPNN策略對模型殘差進行建模與預報處理[24]；

方案12：基于方案6，采用經(jīng)驗半變異函數(shù)模型計算的模型殘差協(xié)方差陣對式(4)中的權(quán)陣非對角線元素進行更新處理，同時采用PLS+BPNN策略對模型殘差進行建模與預報處理。

為對本文提出的改進的鐘差預報模型進行全面分析，將實驗中12個預報方案結(jié)果與WHU的快速鐘差產(chǎn)品進行對比。為討論不同方案的鐘差預報精度，圖5中選取了6套方案中C14與C24衛(wèi)星對比分析建模殘差的差異，通過模型殘差可明顯發(fā)現(xiàn)方案11的建模效果最優(yōu)。

圖5 基于2019年第41天的C14和C24不同方案模型殘差Fig.5 Model residuals of different schemes based on DOY 41, 2019 for C14 and C24 satellites

同時，為具體說明方案7～12中不同實驗的鐘差預報精度，表2給出了不同方案下鐘差預報精度平均值。由于超快速鐘差產(chǎn)品以3h時延和6h間隔進行更新發(fā)布[16]，因此，鐘差預報結(jié)果中僅分析了1天內(nèi)18h的鐘差預報精度。

表2 基于傳統(tǒng)方法、星間相關(guān)性和Lasso算法的18h超快鐘差預報精度(ns)及不同方案提升率

基于圖5與表2中各方案的統(tǒng)計結(jié)果，可以得出：1)與方案2相比，考慮BDS-2/BDS-3衛(wèi)星鐘差相關(guān)性對模型殘差的影響不明顯；2)基于廣泛使用的周期項選取方法進行鐘差建模，將導致BDS-2與BDS-3預報鐘差精度分別下降4.3%和21.1%；3)通過灰色模型對模型殘差進行處理，BDS-2與BDS-3鐘差預報精度分別提高了23.3%與16.9%。需要說明的是，方案8是在傳統(tǒng)方法的基礎(chǔ)上進行建模優(yōu)化處理，需對其模型殘差進行補償。因此，實驗中設(shè)置了基于BPNN和PLS+BPNN算法的鐘差殘差序列預報方案。實驗結(jié)果表明，PLS+BPNN策略可實現(xiàn)模型殘差較傳統(tǒng)BPNN算法建模降低0.5%和6.2%。同時，由于模型殘差序列中存在時空相關(guān)性，在系數(shù)估計中引入半變異函數(shù)對權(quán)陣進行精化處理，方案12較方案11可進一步實現(xiàn)BDS-2與BDS-3預報鐘差精度分別提升8.0%與11.1%。

3 結(jié)論

本文基于前期研究，對BDS-2/BDS-3聯(lián)合超快速鐘差預報策略進行了優(yōu)化處理，并通過3組共12套實驗方案進行驗證，可以得出：

1)由于傳統(tǒng)鐘差模型的兩步求解過程存在精度損失現(xiàn)象，提出了一步估計模型參數(shù)(趨勢項和周期項)的策略；并引入機器學習中的稀疏建模方法對鐘差模型參數(shù)進行自動篩選；基于鐘差預報結(jié)果可知，一步建模策略可略微提升鐘差預報精度。

2)考慮到BDS衛(wèi)星(BDS-2與BDS-3)的差異及星載原子鐘的特點，文中將星間相關(guān)性作為一種間接方法增強模型系數(shù)求解；結(jié)果表明，基于一步建模策略，通過引入BDS-2與BDS-3星間相關(guān)性，可實現(xiàn)BDS-2與BDS-3衛(wèi)星18h鐘差預報精度分別提升27.2%與28.6%；

3)在建立鐘差預報模型時，基于鐘差殘差序列，構(gòu)造經(jīng)驗半變異函數(shù)模型提取殘差序列中的時間相關(guān)性，并將其納入模型系數(shù)估計的權(quán)陣更新中。與已有文獻提出的方法相比[6,24]，該策略可進一步實現(xiàn)BDS-2與BDS-3衛(wèi)星鐘差預報精度分別提升8.0%和11.1%。

綜上，本文提出的BDS-2/BDS-3聯(lián)合超快速衛(wèi)星鐘差預報策略可改善分析中心北斗鐘差產(chǎn)品的精度。但是，本文改進后的策略只基于1個月的實驗進行測試分析，因此，針對BDS-2/BDS-3聯(lián)合超快速衛(wèi)星鐘差預報策略的可用性和準確性還需進一步研究和分析。