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

1. 中國(guó)科學(xué)院 國(guó)家授時(shí)中心，西安 710600

2. 中國(guó)科學(xué)院 精密導(dǎo)航定位與定時(shí)技術(shù)重點(diǎn)實(shí)驗(yàn)室，西安 710600

3. 中國(guó)科學(xué)院大學(xué)，北京 100049

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

1 引言

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

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

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

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

2 N-BEATS算法

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

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

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

3 ZTD預(yù)報(bào)試驗(yàn)設(shè)計(jì)

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

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

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

3.2 基于N-BEATS的ZTD預(yù)報(bào)策略

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

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

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

圖3 預(yù)報(bào)策略示意Fig.3 Schematic diagram of forecast strategy

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

圖4 訓(xùn)練集與驗(yàn)證集示意Fig.4 Schematic diagram of training set and validation set

4 ZTD預(yù)報(bào)精度評(píng)估

4.1 重復(fù)性評(píng)估

為了衡量該算法的魯棒性，針對(duì)BJFS站ZTD，使用3種策略，各重復(fù)訓(xùn)練了9組模型，對(duì)該算法進(jìn)行了重復(fù)性評(píng)估。

選取了連續(xù)365 d的ZTD預(yù)報(bào)值進(jìn)行預(yù)報(bào)精度的重復(fù)性評(píng)估。圖5為該365 d預(yù)報(bào)序列的示意。弧段從2018年7月1日至2019年6月30日。以24 h為邊界，將每天0點(diǎn)的24 h預(yù)報(bào)值拼接為365 d的預(yù)報(bào)值。

圖5 365 d預(yù)報(bào)序列示意Fig.5 Schematic diagram of 365 days forecast

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

表2 9次365d預(yù)報(bào)均值及其STD的統(tǒng)計(jì)結(jié)果

圖 6 三種預(yù)報(bào)策略的均值與STDFig.6 Mean and standard deviation of 3 forecast strategies

4.2 不同弧段ZTD預(yù)報(bào)精度評(píng)估

第4.1小節(jié)的結(jié)果表明，該算法具有較為穩(wěn)健的重復(fù)性。本小節(jié)針對(duì)不同策略，對(duì)選取的9個(gè)IGS跟蹤站的ZTD，進(jìn)行不同預(yù)報(bào)弧段的精度評(píng)估。由于預(yù)報(bào)策略為每2 h進(jìn)行一組預(yù)報(bào)，且每組以2 h為間隔，預(yù)報(bào)24 h，將每組的第0、2、4、…、22小時(shí)的2 h預(yù)報(bào)數(shù)據(jù)分別進(jìn)行拼接得到2018年7月1日至2019年6月30日共365 d的預(yù)報(bào)值。示意如圖7所示，圖中以第2小時(shí)的2 h預(yù)報(bào)值所拼接的365 d預(yù)報(bào)時(shí)序?yàn)槔?，稱之為4 h 365 d預(yù)報(bào)時(shí)序。

圖7 4 h 365 d預(yù)報(bào)序列示意Fig.7 Schematic diagram of 365 days forecast of the 4th hour

(1)殘差序列均值

對(duì)9個(gè)跟蹤站的ZTD使用S1、S2、S3三種策略的不同弧段拼接所得的365 d預(yù)報(bào)值的預(yù)報(bào)殘差均值進(jìn)行評(píng)估。3種策略不同弧段的預(yù)報(bào)殘差均值的統(tǒng)計(jì)結(jié)果如圖8所示。

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

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

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

表3 九個(gè)站不同預(yù)報(bào)弧段預(yù)報(bào)殘差絕對(duì)值的均值

(2)殘差序列的標(biāo)準(zhǔn)差

本小節(jié)對(duì)9個(gè)跟蹤站3種策略的ZTD預(yù)報(bào)殘差序列的標(biāo)準(zhǔn)差進(jìn)行評(píng)估。圖9為策略S1、S2、S3對(duì)應(yīng)的不同預(yù)報(bào)弧段的預(yù)報(bào)殘差STD統(tǒng)計(jì)結(jié)果。

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

從圖9中可以看出，隨著預(yù)報(bào)弧長(zhǎng)的增加，ZTD預(yù)報(bào)殘差的STD隨之增加；但隨著用于預(yù)報(bào)的輸入弧長(zhǎng)增加，預(yù)報(bào)殘差的STD并沒有明顯的改善趨勢(shì)。表4將所有站不同預(yù)報(bào)弧段預(yù)報(bào)殘差STD的均值進(jìn)行了統(tǒng)計(jì)?？筛逦乜吹?，24 h預(yù)報(bào)殘差STD的均值小于30 mm，6 h的預(yù)報(bào)殘差STD的均值小于13 mm。但3種策略不同弧段的預(yù)報(bào)殘差STD均在同一量級(jí)。

表4 九個(gè)站ZTD不同預(yù)報(bào)弧段預(yù)報(bào)殘差STD的均值

5 分析與討論

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

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

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

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

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

6 結(jié)論

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

首先對(duì)N-BEATS算法的重復(fù)性進(jìn)行了評(píng)估。從BJFS站3種策略的9組預(yù)報(bào)殘差結(jié)果來(lái)看，該算法具有良好的重復(fù)性。

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

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

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