基于N-BEATS的單站對流層天頂總延遲預報

蘇行，楊韜，孫保琪，楊旭海

1. 中國科學院 國家授時中心，西安 710600

2. 中國科學院 精密導航定位與定時技術重點實驗室，西安 710600

3. 中國科學院大學，北京 100049

4. 西北工業大學 無人系統技術研究院，西安 710072

1 引言

全球導航衛星系統(global navigation satellite system，GNSS)參數估計中，對流層延遲與測站鐘差、測站高程強相關[1]，是GNSS數據處理中最重要的誤差源之一。信號傳輸過程中，由對流層引起的信號時延一般稱為對流層延遲。不同高度角的對流層延遲通常可由對流層天頂延遲與相應的映射函數計算而來。對流層天頂總延遲(zenith total delay，ZTD)主要由天頂干延遲(zenith hydrostatic delay，ZHD)和天頂濕延遲(zenith wet delay，ZWD)兩部分組成[2]。與水汽強相關的ZWD變化幅度較大，而ZHD較為穩定。因此在實際應用中，通常將ZWD作為未知參數來進行估計[3]。而高精度的對流層延遲先驗值有助于參數估計的快速收斂[4]，從而提升其在實時應用中的可能性[5]。

目前，通過使用事后精密產品進行精密單點定位解算，國際GNSS服務(international GNSS service，IGS)可提供數據間隔為300 s的對流層ZTD事后精密產品。包含350余個IGS跟蹤站的對流層ZTD及其水平梯度分量，產品延遲約為3周[6]。維也納科技大學基于高精度ECMWF氣象數據庫[7]，提供數據間隔為6 h，以地表為基準的VMF1[8]及VMF3[9]事后格網對流層產品，包含ZHD、ZWD及其映射函數系數，其中VMF1的分辨率為2°×2.5°，VMF3的水平分辨率為1°×1°，該格網產品的延遲為1 d。

為了滿足實時、近實時的需求，通常使用經驗模型來進行先驗對流層延遲建模，常見的有GPT及其映射函數GMF[10]、GPT2[11]、GPT2w[12]、UNB3m[13]等。此外，基于ECMWF高精度預報產品，維也納科技大學提供了對流層預報格網產品VMF1-FC[14]，預報范圍為42 h，數據間隔為6 h。

大氣中水汽含量變化迅速且幅值不定，因此對流層產品的時間分辨率越高，越適合于實時應用。利用傳統回歸方法可設計手工特征學習有效模式，而基于數據驅動的深度學習，通過大量有標簽樣本，可端到端進行模式學習，規避傳統建模中的手工特征提取，即基于大量訓練樣本，學習較難建模的特征。本文基于高精度高分辨率氣象數據庫，利用深度學習N-BEATS算法，進行時間分辨率為2 h的對流層ZTD短期預報研究，并進行了精度評估。

2 N-BEATS算法

經典的時間序列預測方法，例如自回歸滑動平均模型(ARMA)和差分整合滑動平均自回歸模型(ARIMA)，要求時間序列數據具有平穩性或者微分后是平穩的。 本質上，此類方法只捕獲線性關系。近年來，隨著算法的進步和計算機算力的提升，深度學習方法強大的自動模式學習能力和非線性函數擬合能力逐漸凸顯，在許多領域的性能都優于傳統方法，例如計算機視覺、語音識別、自然語言處理和時間序列預測。包括卷積神經網絡(CNN)[15-17]、遞歸神經網絡(RNN)[18]和注意力機制[19]在內的經典神經網絡模塊均被用于時間序列預測。

最近，基于神經網絡基底擴展分析的可解釋的時間序列預測(neural basis expansion analysis for interpretable time series forecasting，N-BEATS)在幾個著名的公開數據集上顯示出優于其他時間序列預測方法的性能[20]。N-BEATS的網絡結構可以快速訓練并使神經網絡具有一定的可解釋性。

圖1 N-BEATS網絡結構Fig.1 The network architecture of N-BEATS

3 ZTD預報試驗設計

3.1 數據源

本文使用的高精度高分辨率對流層延遲數據來源于捷克Pecny天文臺GOP-TropDB對流層氣象數據庫[21]。GOP-TropDB基于ECMWF的再分析數據庫ERA-Interim[7]，為用戶提供對流層相關氣象參數，避免了用戶直接從ECMWF獲取氣象參數時面臨的巨大數據量及復雜的數據處理。由于深度學習的模型訓練特點是基于大量數據來進行特征提取，數據量越大越有利于特征的提取。因此本文選取了較長數據弧段，即2002年1月1日至2019年6月30日，采樣間隔為2 h，研究對象為9個在全球分布較為均勻的IGS跟蹤站，如圖2所示。其中，HARB和KOKB的測站高程高于1 000 m，分別為1 558 m、1 167 m，WTZR高程約為666 m，其余6個測站均低于200 m。

圖2 9個IGS跟蹤站橢球高程及分布Fig.2 Global distribution and ellipsoidal heights of 9 IGS tracking stations

3.2 基于N-BEATS的ZTD預報策略

在模型搭建中，根據實測數據的規律，可基于N-BEATS進行趨勢項及周期項模型的混合搭建。首先針對所選測站的ZTD序列進行快速傅里葉變化，以進行周期項的探測。結果表明，除了周年項、半周年項此類長周期項以外，還存在多個短周期項，例如與地球自轉相關的24 h周期項、與月球潮汐相關的12 h短周期項。在本文的模型搭建中，根據不同周期項對各站影響的顯著性，為各站設定不同的周期項個數。此外，經驗證，趨勢項采用默認值即可。

模型搭建好之后，進行輸入及預報步長的設定，本文針對分辨率為2 h的24 hZTD預報，根據輸入弧長設計了3種ZTD預報策略，如表1所示，以2 h為滑動窗口，分別以24 h、48 h、72 h為輸入步長，進行以24 h為步長的預報，示意如圖3所示。

表1 三種策略的輸入及預報弧段設置

圖3 預報策略示意Fig.3 Schematic diagram of forecast strategy

隨后確定學習率，學習率決定了目標模型的收斂速度及精度，經試驗，本文選取的學習率為0.000 1。最后通過大量數據樣本進行模型訓練。以2018年7月1日為分界，將數據分為訓練數據集和驗證數據集。模型的訓練將在訓練集上進行，由此訓練出的模型在驗證集上進行驗證。以圖4中 BJFS站對流層ZTD時間序列為例，取最后一年為驗證集(紅線)，時段為2018年7月1日至2019年6月30日，其余時段(藍線)為訓練集。

圖4 訓練集與驗證集示意Fig.4 Schematic diagram of training set and validation set

4 ZTD預報精度評估

4.1 重復性評估

為了衡量該算法的魯棒性，針對BJFS站ZTD，使用3種策略，各重復訓練了9組模型，對該算法進行了重復性評估。

選取了連續365 d的ZTD預報值進行預報精度的重復性評估。圖5為該365 d預報序列的示意。弧段從2018年7月1日至2019年6月30日。以24 h為邊界，將每天0點的24 h預報值拼接為365 d的預報值。

圖5 365 d預報序列示意Fig.5 Schematic diagram of 365 days forecast

將ZTD預報值與真值相減得到殘差序列，表2為S1、S2、S3三種策略各自所訓練的9組模型的365d殘差序列的均值及其STD(standard deviation)。圖6為3種策略各9組殘差序列均值(左)和STD(右)的箱圖。可以看出，殘差序列的均值大多在亞毫米量級，STD大多約為27～28 mm。由此可見，該算法具有較為穩定的重復性。

表2 9次365d預報均值及其STD的統計結果

圖 6 三種預報策略的均值與STDFig.6 Mean and standard deviation of 3 forecast strategies

4.2 不同弧段ZTD預報精度評估

第4.1小節的結果表明，該算法具有較為穩健的重復性。本小節針對不同策略，對選取的9個IGS跟蹤站的ZTD，進行不同預報弧段的精度評估。由于預報策略為每2 h進行一組預報，且每組以2 h為間隔，預報24 h，將每組的第0、2、4、…、22小時的2 h預報數據分別進行拼接得到2018年7月1日至2019年6月30日共365 d的預報值。示意如圖7所示，圖中以第2小時的2 h預報值所拼接的365 d預報時序為例，稱之為4 h 365 d預報時序。

圖7 4 h 365 d預報序列示意Fig.7 Schematic diagram of 365 days forecast of the 4th hour

(1)殘差序列均值

對9個跟蹤站的ZTD使用S1、S2、S3三種策略的不同弧段拼接所得的365 d預報值的預報殘差均值進行評估。3種策略不同弧段的預報殘差均值的統計結果如圖8所示。

從圖8中可以看出，策略S1中，有個別跟蹤站諸如BJFS、HARB、KOUR的殘差均值隨著預報弧段的增加沒有明顯的起伏變化，但整體上隨著預報弧段的增加，9個跟蹤站的殘差均值呈發散趨勢，6 h以內的預報殘差均值可小于1 mm。策略S2中，除了YAR2站的預報殘差均值隨著預報弧段的增加有明顯的發散趨勢，其余跟蹤站的預報殘差均值的波動大約在2 mm以內。策略S3中，隨著預報弧段增加，YAR2站的預報殘差均值依然比較發散，其余跟蹤站預報殘差均值的波動大多在1 mm以內。

圖8 分別使用S1、S2以及S3策略所得的第N小時預報值的全年殘差均值Fig.8 Mean of the Nth hour forecast residuals of 365 days using strategy S1, S2 and S3

由此可見，隨著預報弧段的增加，預報殘差的均值越發散；用于預報的輸入弧段越短，隨著預報弧長的增加，預報弧段殘差均值越發散。策略S3性能良好，但同時需要較多數據量。表3為所有站不同預報弧段的預報殘差絕對值的均值，可以看出，總體上，3種策略的24 h內預報值均小于3 mm，12 h內的預報值大多優于1 mm；6 h內的預報精度均可達亞毫米量級。

表3 九個站不同預報弧段預報殘差絕對值的均值

(2)殘差序列的標準差

本小節對9個跟蹤站3種策略的ZTD預報殘差序列的標準差進行評估。圖9為策略S1、S2、S3對應的不同預報弧段的預報殘差STD統計結果。

圖9 分別使用S1、S2以及S3策略所得的第N小時預報值的全年殘差STDFig.9 STD of the Nth hour forecast residuals of 365 days using strategy S1, S2 and S3

從圖9中可以看出，隨著預報弧長的增加，ZTD預報殘差的STD隨之增加；但隨著用于預報的輸入弧長增加，預報殘差的STD并沒有明顯的改善趨勢。表4將所有站不同預報弧段預報殘差STD的均值進行了統計。可更清晰地看到，24 h預報殘差STD的均值小于30 mm，6 h的預報殘差STD的均值小于13 mm。但3種策略不同弧段的預報殘差STD均在同一量級。

表4 九個站ZTD不同預報弧段預報殘差STD的均值

5 分析與討論

由于對流層延遲變化頻繁，為了更直觀地進行分析，對9個站的對流層ZTD、ZHD及ZWD進行了二次多項式擬合并繪圖，擬合弧段從2018年7月1日至2019年6月30日共365 d，如圖10所示。

從圖10中可以看出，ZHD較為平穩且幅值較大，決定了ZTD的幅值量級，ZWD幅值比ZHD小但變化幅度較大，具有較為明顯的季節性趨勢，是構成ZTD波動的主要因素。除HARB和KOKB站以外，其余的ZHD幅值均約為2.2 m，結合9個站的地理位置分布(見圖2)，可以看出測站高程與ZHD幅值成反比。位于赤道附近的KOUR站的ZTD幅值雖大，但波動范圍適中，其預報殘差序列的均值及STD均處于中間水平。位于南極洲的OHI2站、MCM4站的ZTD幅值及波動較小，其預報殘差的均值及STD均處于較好水平。可見，預報精度與ZWD的波動密切相關。

圖10 2018年7月至2019年6月9個IGS跟蹤站的對流層的平滑曲線ZTD，ZHD，ZWDFig.10 Smoothed tropospheric ZTD, ZHD and ZWD of 9 IGS tracking stations from Jul. 2018 to Jun. 2019

隨著輸入步長的增加，預報殘差的均值均呈降低趨勢。但3種策略的預報殘差STD均在同一量級。經觀察，該量級與ZWD的變化范圍有關。對使用策略S1的9個站的不同弧段ZTD預報殘差的STD求均值，然后對9個站的STD均值和其各自對應的ZWD變化范圍求相關性系數。如圖11所示，預報殘差的STD與ZWD變化范圍的相關性系數為0.85。作為影響ZTD波動的主要因素，ZWD變化頻繁、較難捕捉規律，導致其仍是限制ZTD預報精度的主要原因。而預報殘差所對應的是數據中未能進行有效建模的部分，因此預報殘差STD與ZWD變化范圍具有強相關。結合圖2中的地理分布及圖10中的ZWD幅值，可以看出ZWD幅值與緯度強相關，并且赤道與極地范圍內的ZWD變化范圍較小。因此，在后續的模型改進中，可考慮將緯度作為重要因素引入模型。

圖11 九個站ZWD與預報殘差STD的相關性示意Fig.11 Correlation between ZWD and standard deviation of forecast residuals of 9 stations

6 結論

基于N-BEATS數據驅動算法，設計了3種對流層ZTD預報策略，用于預報的輸入數據弧長分別為24 h、48 h和72 h。選取了9個分布均勻的IGS跟蹤站進行研究，數據弧段從2002年1月1日到2019年6月30日共計18.5 a，選取前17.5 a來訓練模型，最后一年用來進行精度評估驗證。

首先對N-BEATS算法的重復性進行了評估。從BJFS站3種策略的9組預報殘差結果來看，該算法具有良好的重復性。

然后對N-BEATS算法的不同弧段預報精度進行了評估。統計結果表明，ZTD的波動范圍決定預報殘差的均值及STD。12 h以內預報殘差的均值大多在亞毫米量級。隨著預報弧長的增加，殘差均值的精度逐漸衰減。但隨著用于預報的輸入弧長的增加，殘差均值的精度衰減可得到有效改善。預報殘差STD的幅值隨著預報弧長的增加，呈增大趨勢，但不同策略的預報殘差STD沒有明顯區別。經觀察分析，預報殘差的STD與ZWD變化范圍有較強相關性，而ZWD變化范圍與緯度存在相關性，因此，在未來的預報研究中將側重于通過引入緯度因子來改善ZWD的影響。

基于N-BEATS的預報策略的12 h以內的平均預報精度可達亞毫米量級，且具有較高時間分辨率，有利于進行諸如實時精密單點定位等GNSS實時應用。

致謝：感謝中國科學院國家授時中心iGMAS 分析中心、國家科技基礎條件平臺-國家空間科學數據中心(http:∥www.nssdc.ac.cn)。