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

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

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

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

4. 西北工業(yè)大學 無人系統(tǒng)技術(shù)研究院，西安 710072

1 引言

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

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

為了滿足實時、近實時的需求，通常使用經(jīng)驗?zāi)Ｐ蛠磉M行先驗對流層延遲建模，常見的有GPT及其映射函數(shù)GMF[10]、GPT2[11]、GPT2w[12]、UNB3m[13]等。此外，基于ECMWF高精度預報產(chǎn)品，維也納科技大學提供了對流層預報格網(wǎng)產(chǎn)品VMF1-FC[14]，預報范圍為42 h，數(shù)據(jù)間隔為6 h。

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

2 N-BEATS算法

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

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

圖1 N-BEATS網(wǎng)絡(luò)結(jié)構(gòu)Fig.1 The network architecture of N-BEATS

3 ZTD預報試驗設(shè)計

3.1 數(shù)據(jù)源

本文使用的高精度高分辨率對流層延遲數(shù)據(jù)來源于捷克Pecny天文臺GOP-TropDB對流層氣象數(shù)據(jù)庫[21]。GOP-TropDB基于ECMWF的再分析數(shù)據(jù)庫ERA-Interim[7]，為用戶提供對流層相關(guān)氣象參數(shù)，避免了用戶直接從ECMWF獲取氣象參數(shù)時面臨的巨大數(shù)據(jù)量及復雜的數(shù)據(jù)處理。由于深度學習的模型訓練特點是基于大量數(shù)據(jù)來進行特征提取，數(shù)據(jù)量越大越有利于特征的提取。因此本文選取了較長數(shù)據(jù)弧段，即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預報策略

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

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

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

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

隨后確定學習率，學習率決定了目標模型的收斂速度及精度，經(jīng)試驗，本文選取的學習率為0.000 1。最后通過大量數(shù)據(jù)樣本進行模型訓練。以2018年7月1日為分界，將數(shù)據(jù)分為訓練數(shù)據(jù)集和驗證數(shù)據(jù)集。模型的訓練將在訓練集上進行，由此訓練出的模型在驗證集上進行驗證。以圖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組模型，對該算法進行了重復性評估。

選取了連續(xù)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。由此可見，該算法具有較為穩(wěn)定的重復性。

表2 9次365d預報均值及其STD的統(tǒng)計結(jié)果

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

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

第4.1小節(jié)的結(jié)果表明，該算法具有較為穩(wěn)健的重復性。本小節(jié)針對不同策略，對選取的9個IGS跟蹤站的ZTD，進行不同預報弧段的精度評估。由于預報策略為每2 h進行一組預報，且每組以2 h為間隔，預報24 h，將每組的第0、2、4、…、22小時的2 h預報數(shù)據(jù)分別進行拼接得到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種策略不同弧段的預報殘差均值的統(tǒng)計結(jié)果如圖8所示。

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

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

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

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

(2)殘差序列的標準差

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

圖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的均值進行了統(tǒng)計?？筛逦乜吹?，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較為平穩(wěn)且幅值較大，決定了ZTD的幅值量級，ZWD幅值比ZHD小但變化幅度較大，具有較為明顯的季節(jié)性趨勢，是構(gòu)成ZTD波動的主要因素。除HARB和KOKB站以外，其余的ZHD幅值均約為2.2 m，結(jié)合9個站的地理位置分布(見圖2)，可以看出測站高程與ZHD幅值成反比。位于赤道附近的KOUR站的ZTD幅值雖大，但波動范圍適中，其預報殘差序列的均值及STD均處于中間水平。位于南極洲的OHI2站、MCM4站的ZTD幅值及波動較小，其預報殘差的均值及STD均處于較好水平?？梢?，預報精度與ZWD的波動密切相關(guān)。

圖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均在同一量級。經(jīng)觀察，該量級與ZWD的變化范圍有關(guān)。對使用策略S1的9個站的不同弧段ZTD預報殘差的STD求均值，然后對9個站的STD均值和其各自對應(yīng)的ZWD變化范圍求相關(guān)性系數(shù)。如圖11所示，預報殘差的STD與ZWD變化范圍的相關(guān)性系數(shù)為0.85。作為影響ZTD波動的主要因素，ZWD變化頻繁、較難捕捉規(guī)律，導致其仍是限制ZTD預報精度的主要原因。而預報殘差所對應(yīng)的是數(shù)據(jù)中未能進行有效建模的部分，因此預報殘差STD與ZWD變化范圍具有強相關(guān)。結(jié)合圖2中的地理分布及圖10中的ZWD幅值，可以看出ZWD幅值與緯度強相關(guān)，并且赤道與極地范圍內(nèi)的ZWD變化范圍較小。因此，在后續(xù)的模型改進中，可考慮將緯度作為重要因素引入模型。

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

6 結(jié)論

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

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

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

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

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