王迪,王茂磊,楊宇飛,李海航,葛條
( 1. 北京衛星導航中心, 北京 100089;2. 北京開運聯合信息技術集團股份有限公司, 北京 100020;3. 北京工業大學, 北京 100124 )
GNSS 是指能在地球表面或近地空間的任何地點為用戶提供全天候三維坐標和速度以及時間信息的空基無線電導航定位系統[1].目前全球有四大衛星導航系統:美國的GPS、俄羅斯的GLONASS、歐盟的Galileo 和中國的北斗衛星導航系統(BeiDou Navigation Satellite System,BDS).BDS 經過“三步走”發展戰略,實現服務性能大幅提升.根據BDS 公開服務性能規范(3.0 版)定義,定位、導航和授時(positioning,navigation and timing,PNT)服務精度包括定位精度、測速精度、授時精度和服務可用性等[2].
近年來,機器學習技術在多個領域得到了廣泛應用,在GNSS 方面也有一定的研究進展.Hsu[3]采用機器學習的方法探討了多路徑效應檢測;Linty 等[4]根據機器學習決策樹算法研究了GNSS 中的電離層閃爍;周相兵[5]利用城市出租車GNSS 數據,采用智能聚類學習算法對數據進行了分類,為城市道路規劃和路網更新提供了新方法;Xia 等[6]研究了基于混合機器學習系統的車輛GNSS 觀測異常檢測.上述研究均為利用機器學習實現了對GNSS 數據進行分類或聚類.此外,也有學者基于機器學習,對觀測結果或某種現象進行預測.如Kiani[7]將機器學習方法應用于GNSS 時間序列研究,對地面沉降或隆起進行了預測;駱黎明[8]基于樹模型反演了GNSS-R 海面風場.盡管機器學習在上述方面的應用百花齊放,但基于機器學習的無源定位導航服務性能評估卻鮮有人研究.
眾所周知,傳統數理方法的GNSS 性能評估是基于大地測量與導航測繪原理開展,其所有數據依賴于地面測站的觀測與傳輸,導致地面測站的布局對定位精度的評估、預測存在一定的關聯性和局限性,而機器學習可以較好地挖掘這方面的信息,通過對大量數據進行篩選、特征提取、模型訓練、歸納總結等,探索研究機器學習在導航定位服務中的應用效果.
本文利用機器學習算法,結合大地測量和導航測繪原理,以與定位精度有關的數據作為特征,對定位精度進行擬合,為定位性能的評估與預測提出了一種新方法和思路.即運用特征提取工程模塊自動提取并合理篩選出影響定位精度的相關數據特征,然后將上述特征輸入機器學習工程模塊進行模型訓練和性能測試,通過特征數據擬合定位精度真值,實現定位精度預測,使用定位精度預測值對服務性能進行大致評估.同時,本文工作本身更多的是基于“傳統統計方法”獲得的定位精度結果進行模型訓練的,則傳統方法中存在的各種“關聯性和局限性”問題在建立的模型中只會放大而不會消失.因此,本文更多的是提供一種新的研究思路,即利用機器學習的思想結合GNSS 理論實現對定位精度的評估和預測.
本文選用武漢大學國際GNSS 服務(International GNSS service,IGS)網站2020—2021 兩年全球測站的觀測文件和每日廣播星歷作為原始數據.選用BDS B1I 頻段信號,以每30s 的時間間隔對測站位置進行單點定位,并與測站位置參考值進行對比,獲得東(east,E)、北(north,N)、天頂(up,U)三方向上的定位誤差σE、σN、σU.
本文將測站實際定位精度作為標簽值,用于給機器學習模型得出的預測值提供參考.
定位精度是一段時間內定位誤差的統計值.在2.1 節中描述了計算E、N、U 三個方向上的定位誤差的方法.接著,經過下面兩步的處理將三個方向上的定位誤差轉換為水平、垂直方向定位精度.
首先將三個方向的誤差轉換為水平、垂直兩個方向上的定位誤差,轉換方法如下:
式中:σH、σV為水平方向誤差與垂直方向誤差;σE、σN、σU為E、N、U 三個方向上的誤差.
其次,將得到的垂直方向、水平方向定位誤差按一定采樣時長進行分段,使用均方差公式對每段時間內的定位誤差進行統計,得到的統計值作為每日不同時間段內的垂直方向、水平方向定位精度.
在本文中,主要采用1h 和3h 的采樣時長進行分段統計.對于采樣時長為1h 的數據集,一天內共有24 組統計值,代表24 個時間段內的定位精度.對于采樣時長為3h 的數據集,一天內共有8 組統計值,代表8 個時間段內的定位精度.
考慮到定位精度可能受到天氣、地理位置等因素的影響,將衛星幾何構型、信號傳播段和信號用戶段等對應數據作為重要特征值,將其按時間順序添加到數據集中,與標簽值一一對應.
衛星信號在傳播時會受到傳播路徑上多種因素的干擾,傳播段上的主要影響因素為電離層狀態和對流層狀態:對于電離層總體狀態,主要受太陽活動與地球地磁活動的影響,因此考慮可反應地球地磁活動情況和太陽活動情況的數據作為特征值;局部的電離層狀態則與當地太陽光照情況有關,可由當地時間間接反應,因此考慮將定位時間作為特征值.在本文中,地球地磁活動情況和太陽活動情況主要通過德國地磁中心提供的地磁指數文件獲得,選取文件中的地磁指數Ap、太陽黑子數SN、太陽射電輻射通量F10.7作為特征值.時間特征值則通過將各個標簽值對應的時間分解為定位時間的年積日(day of year,DOY)以及當日小時數(hour of day,HOD)來獲取.此外,通過對國際全球連續監測評估系統(international GNSS Monitoring&Assessment System,iGMAS)網站提供的電離層文件在時間與空間上線性插值,可以直接得到不同定位時間下測站上方的電離層電子總含量(ionospheric total electron content,TEC)值.上述選取特征都直接或間接地反應了電離層狀態;對流層狀態與氣象情況相關,如定位地點周邊溫度T、大氣壓Po、相對濕度U等,由于氣象參數的獲取受到氣象站分布的限制,許多測站無法獲取周邊準確氣象參數,故不作為特征值考慮;而用戶段的影響包括接收機鐘差、接收機噪聲等,由于不同接收機型號對定位精度評估受人為技術影響較大,故此處不作為特征值考慮.
此外,定位時的衛星幾何構型也將影響定位精度,使用定位時的可見衛星位置計算三維位置精度因子(position dilution of precision,PDOP),并作為一種特征值,可以通過數值大小定量反應衛星的幾何構型的情況.
綜合考量,將衛星幾何構型和信號傳播段的主要因素作為本文機器學習的特征值,具體為PDOP、測站上方的TEC、地球地磁指數Ap、每日太陽黑子數SN、太陽輻射通量F10.7、定位時間DOY 與HOD.
基于GNSS 數據處理具有數據計算量大、數據關系耦合多、計算類別復雜等特點,采取將數據集高維度原始空間投影到低維度特征空間方式,保持樣本類別區分性,降低計算量,減小參數估計誤差從而避免過擬合問題.在獲取特征數據與定位精度數據后,需要將該測站數據按時間拼接起來,構成可供機器學習模型使用的數據集.
將數據輸入機器學習模型之前,需要對數據做預處理,以獲取最佳性能.本文采用周期性數據編碼和數據標準化的處理方式.周期性數據編碼主要使用sin、cos 函數對HOD、DOY 兩種具有周期性的特征進行處理,使得模型可以學習到其周期性.數據標準化則使用以下公式:
式中:為標準化后的第i類數據中的第j個數據;Xij為第i類數據中的第j個數據;μi為第i類數據的均值;σi為第i類數據的標準差.
通過計算不同數據的均值與方差,并使用上式依次進行處理,可以根據不同數據的分布,將所有輸入輸出數據的分布轉化為正態分布.這有助于加速模型收斂速度并提高模型預測準確率.
梯度提升決策樹(gradient boosting decision tree,GBDT)[9]、支持向量回歸(support vector regression,SVR)[10]和多層感知機(multilayer perceptron,MLP)[11]均具有較強非線性擬合能力,且較為常見.其中GBDT是一種具有很強泛化性能的機器學習模型,具有非線性擬合能力,并適用于所有規模的數據集,在多個領域得到廣泛的應用,本文優選GBDT 模型進行定位精度擬合.
GBDT 通過構建多個弱學習器組合得到有較強預測性能的強學習器,其中,弱學習器指預測效果較差的模型,GBDT 采用的弱學習器為回歸決策樹.圖1展示了一個用于回歸任務的GBDT 示意圖,其由M棵回歸決策樹組成.

圖1 GBDT 示意圖
對于圖1 中展示的GBDT 模型,其輸出y的表達式為
式中,βm和Tm(x)分別為第M棵決策樹的輸出權重與輸出值.
GBDT 在訓練時,通過多次迭代逐一生成回歸決策樹并加入到模型中,不斷提高模型預測能力.對于一個由M棵回歸決策樹組成的GBDT 模型來說,需要進行m次迭代生成回歸決策樹.在其中的第m次迭代時,由前m-1 步得到的m-1 棵決策樹組成的強學習器fm-1(x)為
式中,βi和θi分別為第i棵決策樹的權重與參數.此強學習器的預測結果與訓練集中的標簽值的偏差用損失函數L(y,fm-1(x))表示.而對于本次迭代需要生成的決策樹T(x;θm),需要將損失函數L(y,fm-1(x))的負梯度γm作為擬合目標
對于GBDT 回歸模型來說,使用的損失函數L(y,fm(x))為均方差函數
則負梯度rm可以進一步化為
即第m次迭代擬合的決策樹需要以上一個強學習器的預測值與標簽值的差值作為擬合目標.
在得到第m次迭代產生的決策樹T(x;θm)后,將其加入強學習器中
通過重復以上迭代過程,不斷生成新的決策樹并加入強學習器中,直到強學習器中的決策樹數量滿足設定的數量,即完成了整個模型的訓練.
GBDT 完整訓練流程可表示為:
1)初始化強學習模型f0(x)為訓練集標簽值的均值.
2)對m=0,1,2,···,M.
a)根據標簽值和上一個強學習器的預測值計算每個訓練集樣本殘差rmi=yi-fm-1(xi),i=1,2,···,N;
b)用殘差訓練回歸樹,得到T(x;θm);
c)更新當前強學習器fm(x)=fm-1+βmT(x;θm);
3)得到最終的模型fM(x).
將上述數據集,按一定比例分為沒有交集的訓練集與驗證集.其中訓練集用于GBDT 模型的訓練工作,驗證集則用于驗證在未知數據下模型的預測表現.
在訓練開始前,需要進行預處理.對訓練集使用的預處理方法有:針對時間等周期性數據,使用cos、sin 函數映射;對于非數值型數據進行獨熱編碼;對所有數據進行標準差標準化.
在訓練時,需要進行超參數的確定,確定模型在數據集上表現最佳的超參數.本文使用網格搜索的方法確定GBDT 在訓練集上的最佳超參數.首先對整體性能影響最大的決策樹個數進行搜索,在最佳決策樹個數基礎上,對最大樹深、葉子節點最小樣本數等決策樹參數進行搜索,最后在上述最佳超參數基礎上對劃分時考慮的特征數、下采樣率和學習率進行搜索.確定超參數后,使用最佳超參數創建模型,并使用訓練集進行訓練.
為了驗證GBDT 模型可實施性與評估性,需要對GBDT 模型進行測試.首先使用訓練集參數,對驗證集進行標準化處理.然后將驗證集中的特征部分輸入訓練好的模型中進行預測,獲取模型預測值.此處以ENAO 測站水平精度為例闡述模型預測值生成過程.
1)準備輸入特征:根據時間和跟蹤站的經緯高,分別獲取輸入特征Ap、SN、F10.7、TEC、HOD、DOY和PDOP,演示樣本的輸入特征如表1 所示,其中DOY 與HOD 已經過sin、cos 函數映射處理.并使用式(3),根據每類特征的均值與標準差對進行數據標準化,ENAO 訓練集中每類特征的數據分布如表2 所示.根據式(3),對表1 中的每類特征分別減去每類特征對應的均值并除以標準差,最終得到的標準化后的輸入特征如表3 所示.

表1 演示樣本的輸入特征

表2 ENAO 的訓練集分布

表3 處理后的輸入特征
2)使用模型預測:加載ENAO 的水平定位精度模型,并獲取預測值.模型使用實驗中最佳模型GBDT 訓練得到.對于GBDT 模型,其由若干決策回歸樹構成.GBDT 的預測過程就是將其中所有的決策回歸樹的輸出乘以一定權重后累加的過程.對于演示模型,ENAO 的水平精度預測模型而言,其由400 個決策回歸樹構成,這里挑選其中第1 個、第200 個、第400 個決策回歸樹,來演示GBDT 模型中的詳細預測過程.第1 個、第200 個、第400 個決策回歸樹的結構如圖2~4 所示.

圖2 第1 個決策樹的結構

圖3 第200 個決策樹的結構

圖4 第400 個決策樹的結構
將樣本輸入決策樹后,根據樹中每個結點的條件和樣本的數據,將滿足條件的樣本分至左節點,不滿足條件的分至右節點.不斷重復此過程,直至樣本落入葉子節點.最后以葉子節點上的值value 作為此樣本的預測值輸出.對于上述的三個決策樹,將三個樣本分至葉子節點的過程如圖5~7 所示.圖中紅色、綠色、藍色箭頭的路徑分別代表樣本1、樣本2、樣本3 落入葉子節點的過程.

圖5 演示樣本在第1 個決策樹上的輸出過程

圖6 演示樣本在第200 個決策樹上的輸出過程

圖7 演示樣本在第400 個決策樹上的輸出過程
將每個樹的輸出值乘以權重后進行累加,即可得到GBDT 的最終輸出.GBDT 模型參數中的學習率即為權重.演示模型的學習率為0.1,其輸出即為400 個樹的結果累加后乘以0.1.
假設模型僅由展示的三個決策樹組成,則結果為:
樣本1:0.1×(-0.234+0.013-0.004)=-0.0225
樣本2:0.1×(-0.074+0.013-0.008)=-0.0069
樣本3:0.1×(1.27+0.013-0.008)=1.275
此結果為3 個演示樣本在假設模型上的輸出值.
3)處理輸出值:模型上的輸出值并不能直接作為定位精度的預測值,需要進一步恢復成未標準化時的數值,這個過程稱為反標準化.這是由于模型在訓練時擬合的都是標準化后的標簽值.為了進行反標準化,需要得知數據平均值與標準差,并將數值乘以標準差后再加上平均值.對于演示模型來說,訓練集中水平定位精度平均值為0.902,標準差為0.165.則對于上述三個預測值來說,其反標準化結果為
樣本1:-0.0225×0.165+0.902=0.8982
樣本2:-0.0069×0.165+0.902=0.9008
樣本3:1.275×0.165+0.902=1.1123.
此結果為3 個樣本在假設模型上對水平定位精度的預測值.預測值將從機器學習模塊中返回.
為了評估模型的預測性能,本文主要使用1-MAPE作為模型預測準確率.其中MAPE 代表平均絕對百分比誤差,其計算公式為
式中:ylabel為樣本的標簽值;ypred為樣本的預測值.
MAPE 主要衡量誤差絕對值與真實值之間的比值,反應預測值與標簽值之間的不符合程度,其值域為(0,∞).標簽值與預測值越相符,誤差越小,則MAPE 越接近0;反之,越不相符,誤差越大,MAPE越接近正無窮.
而使用1-MAPE 作為模型預測準確率指標,可以直觀反應預測值與標簽值之間的符合程度,其值域為(-∞,1).標簽值與預測值越接近時,預測準確率1-MAPE 越接近1;反之越接近負無窮.
在評估性能時,首先計算所有測試樣本的MAPE的平均值作為模型的MAPE 指標值.其次計算1-MAPE作為模型的預測準確率.
將測站數據分為訓練模型學習能力的訓練集和評估模型泛化能力的驗證集.訓練集和驗證集的數據劃分原則:在數據集中按一定比例隨機抽取,但需保證訓練集和驗證集兩個數據集合互斥.
實驗環境:計算服務器配置為i9-12900K 和RTX3090 顯卡,利用Python 語言在數據分析處理包Pandas 上進行數據集加載、劃分等,在機器學習包scikit-learn 進行模型創建、數據標準化、模型訓練與測試、模型保存等.
使用ABMF、JFNG、MKEA 等11 個測站數據組成的小型數據集,對數據集上的水平定位精度進行擬合計算與評估,從而比較GBDT、支持向量回歸(support vaetor regression,SVR)和多層感知器(multi-layer perceptron,MLP)在本任務中的性能.其中,由于MLP對數據量要求較大,使用的是采樣時間為1h 的數據集,GBDT 和SVR 使用的是采樣時間為3h 的數據集.將數據集按8∶2 的比例分為訓練集與驗證集,根據驗證集中樣本預測準確率指標的直方分布圖比較模型性能.三種模型的測試結果如圖8~10,其中橫坐標代表預測準確率,縱坐標代表相應樣本數量,上方數字代表驗證集的數據量.

圖8 GBDT 測試結果
如圖8 所示,GBDT 的測試結果如下: 預測準確率處于0.9~1 的測試樣本占總樣本的21.5%,處于0.8~0.9 的測試樣本占20.2%,處于0.7~0.8 的測試樣本占16.9%.如圖9 所示,SVR 處于0.9~1 的測試樣本占總樣本的20.8%,處于0.8~0.9 的測試樣本占18.2%,處于0.7~0.8 的測試樣本占16%.如圖10 所示MLP 處于0.9~1 的測試樣本占總樣本的18.4%,處于0.8~0.9 的測試樣本占16.8%,處于0.7~0.8 的測試樣本占14.3%.

圖9 SVR 測試結果

圖10 MLP 測試結果
預測準確率的值越接近1,代表樣本的預測值與標簽值越接近.因此,預測準確率值越靠近1 的測試樣本在總測試樣本中的占比越大代表預測效果越好.上述結果表明,三種模型對定位精度預測性能的排名為GBDT>SVR>MLP.另外,模型訓練時間排序為SVR>GBDT>MLP.綜合考慮預測性能與訓練時間,結果表明本文所選GBDT 模型最適合定位性能評估任務.
使用GBDT 模型對DGAR、MIZU、JFNG、CUSV等140 余個測站數據集,分別對單測站完成建模,將每個測站的數據集按9∶1 的比例分為訓練集與驗證集,得出訓練結果如下.
在所有測站中,中國及周邊區域12 個測站模型模擬定位精度預測準確率1-MAPE 為92.36%,最差為PTGG 站,預測準確率1-MAPE 為89.26%,全球范圍120 個測站模型模擬定位精度預測準確率1-MAPE 為86.59%,最差為SCOR 站,預測準確率1-MAPE 為81.46%.圖11 為測站評估結果.

圖11 測站GBDT 模型評估結果
結果表明,GBDT 模型用于衛星導航全球定位精度評估效果與傳統數理統計框架下得到的實測值較為吻合,該方法可為后續研究機器學習在基于時間和空間條件下對全球定位性能評估問題提供理論基礎和經驗.
本文提出了一種基于機器學習模型評估衛星導航系統定位性能的方法,主要通過模型對定位精度實現高準確率預測,進而評估定位性能.在相同數據集上,選取三種常用于非線性擬合的機器學習模型進行訓練,得到了GBDT 模型更適合衛星導航定位性能評估的結論.同時,對全球共140 余個測站分別進行了單獨建模,結果表明:機器學習擬合得出的導航定位精度評估效果與實測值較為吻合,說明基于機器學習模型評估衛星導航定位性能的方法可行有效,為下一步對定位性能在時空域的預測提供了技術基礎.
但是,本文采取的方法還存在諸多不足.如特征值只充分考慮了傳播段,對用戶段和空間段考慮不足;模型超參數搜索方法較為簡單等.后續將進一步增加對GNSS 數據相關特征選取方面的研究,以提高評估和預測性能;改進模型超參數搜索方法,使用如遺傳算法(genetic algorithm,GA)等方法尋找模型最佳超參數,避免在超參數搜索時因手動進行網格搜索帶來人為引入的局限性.