一種區(qū)域精密天頂對(duì)流層延遲組合預(yù)報(bào)模型

劉江濤 劉雙童 葉正真 高志鈺

1 甘肅農(nóng)業(yè)職業(yè)技術(shù)學(xué)院，蘭州市段家灘路425號(hào)，730030 2 中國地震局地質(zhì)研究所地震動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京市華嚴(yán)里甲1號(hào)，100029

對(duì)流層延遲指電磁波信號(hào)在通過高度40 km以下未被電離的中性大氣層時(shí)所產(chǎn)生的一種信號(hào)延遲，是GNSS領(lǐng)域主要誤差源之一。對(duì)流層延遲可用天頂對(duì)流層延遲(zenith tropospheric delay, ZTD)和與高度角相關(guān)的映射函數(shù)的乘積表示。GNSS中消除(削弱)對(duì)流層延遲的方法大致可以歸納為隨機(jī)模型、待定參數(shù)估計(jì)、差分法、模型改正等。對(duì)流層延遲改正模型既可以為高精度GNSS測(cè)量提供較精確的對(duì)流層延遲先驗(yàn)值，又可以顯著降低相對(duì)定位差分后的殘余對(duì)流層延遲誤差，因此模型改正成為常用的對(duì)流層延遲處理策略。

早期的對(duì)流層延遲模型主要有Hopfield、Saastamoinen和Black模型，需代入實(shí)測(cè)氣象參數(shù)進(jìn)行計(jì)算，使得模型使用受到限制。鑒于此，部分學(xué)者直接針對(duì)氣象參數(shù)建立時(shí)空模型，以Saastamoinen為基礎(chǔ)衍生出UNB、EGNOS、GPT等一系列無氣象參數(shù)的對(duì)流層延遲模型。隨著全球四維天頂對(duì)流層延遲信息的不斷累積，由回歸參數(shù)直接模擬對(duì)流層延遲時(shí)空變化成為目前主流的對(duì)流層延遲建模方法[1]。TropGrid2模型、SHAO系列模型(SHAO-C、SHAO-G、SHAO-H)、IGGtrop模型、GZTD系列模型均是基于高精度ZTD時(shí)間序列建立的全球數(shù)值模型，相比UNB、EGNOS精度得到進(jìn)一步提高，可達(dá)4 cm。由于全球范圍內(nèi)不同區(qū)域的氣候差異較大，這類全球ZTD模型在局部地區(qū)會(huì)存在明顯偏差[2]。因此，部分學(xué)者依托CORS網(wǎng)基站建立了區(qū)域精密對(duì)流層延遲模型[3-4]，雖然可縮小系統(tǒng)偏差，但精度仍約為4 cm。主要原因在于模型時(shí)間域內(nèi)僅由傅里葉級(jí)數(shù)擬合，只能反映ZTD變化主趨勢(shì)，且較難精細(xì)描述ZTD序列的快速變化特征。另一方面，空間域采用函數(shù)(球諧函數(shù)+指數(shù)函數(shù)、多面函數(shù))擬合或三維網(wǎng)格表征精度有限，且存在參數(shù)過多、模型復(fù)雜、使用不便等問題。此外，通過相關(guān)實(shí)驗(yàn)來驗(yàn)證模型實(shí)用效果的研究較少。

針對(duì)上述問題，本文提出一種結(jié)合小波變換、傅里葉級(jí)數(shù)擬合、AR和SVR的ZTD建模方法。該模型將ZTD時(shí)間序列分解為高、低頻序列分別進(jìn)行建模，可最大程度還原ZTD時(shí)間變化特性；空間部分則采用以統(tǒng)計(jì)學(xué)習(xí)理論為基礎(chǔ)的SVR機(jī)器學(xué)習(xí)算法建立關(guān)于空間位置的非線性模型，且只需確定核函數(shù)類型和3個(gè)參數(shù)即可完成建模。本文基于一定數(shù)量的GNSS基站ZTD數(shù)據(jù)完成建模，并進(jìn)行模型內(nèi)、外符合精度檢驗(yàn)。最后通過偽距單點(diǎn)定位實(shí)驗(yàn)測(cè)試該模型的實(shí)用效果。

1 建模方法及原理

利用區(qū)域內(nèi)GNSS基站觀測(cè)數(shù)據(jù)解算的對(duì)流層延遲序列，可研究該區(qū)域內(nèi)對(duì)流層延遲時(shí)空變化，是建立區(qū)域精密對(duì)流層延遲模型的可靠數(shù)據(jù)源。本文所用數(shù)據(jù)為NOAA(national oceanic and atmospheric administration)CORS網(wǎng)在加州境內(nèi)的部分GNSS基站2016～2018年的RINEX觀測(cè)值文件(https:∥geodesy.noaa.gov/corsdata/rinex/)經(jīng)GAMIT解算得到的ZTD和ZHD(天頂干延遲)序列，其中94個(gè)站點(diǎn)用于建模，24個(gè)站點(diǎn)(靠近海岸線，對(duì)流層延遲變化相對(duì)劇烈)用于模型預(yù)報(bào)測(cè)試。

GAMIT采用雙差觀測(cè)值，需引入長距離測(cè)站來獲取絕對(duì)延遲。本文采用GAMIT進(jìn)行數(shù)據(jù)處理時(shí)引入周邊4個(gè)長距離IGS站(MKEA、BAKE、SASK、WHIT)參與解算。詳細(xì)解算流程參考文獻(xiàn)[5]。

為研究ZTD的時(shí)間變化特征，隨機(jī)選取1個(gè)站點(diǎn)(AZU1)2016～2018年的ZTD時(shí)間序列進(jìn)行快速傅里葉變換，圖1為時(shí)間序列和對(duì)應(yīng)的頻譜圖。

圖1 AZU1站ZTD時(shí)間序列及其頻譜圖Fig.1 ZTD time series and its sequence spectrum of AZU1 station

從圖1(b)可以看出，AZU1測(cè)站ZTD序列含有2個(gè)顯著頻率(0.002 7 Hz和0.005 8 Hz)，其對(duì)應(yīng)周期分別為365.25 d和365.25/2 d，表現(xiàn)出明顯的年周期和半年周期特性，符合ZTD普遍時(shí)間變化規(guī)律。因此大部分學(xué)者在時(shí)間域均采用式(1)，即傅里葉級(jí)數(shù)進(jìn)行最小二乘擬合:

(1)

式中，T1=365.25，T2=365.25/2，doy為年積日，A0(年均值)、A1(年振幅)、A2(半年振幅)、φ1(年相位)、φ2(半年相位)為待擬合參數(shù)。

天頂對(duì)流層延遲包含ZHD和ZWD(天頂濕延遲)，其中ZWD受水汽含量影響變化較為劇烈，導(dǎo)致ZTD整體出現(xiàn)圖1(a)所示的高頻振蕩特性。若采用式(1)直接擬合ZTD序列，擬合殘差較大。文獻(xiàn)[6]利用小波分解技術(shù)提取ZTD序列的低頻部分來研究其變化規(guī)律，結(jié)果表明ZTD低頻序列可反映ZTD序列主變化趨勢(shì)，但文獻(xiàn)未對(duì)高頻序列進(jìn)行分析，而該部分包含ZTD細(xì)節(jié)信息。因此，本文采用AR模型剔除噪聲，進(jìn)一步提取這部分細(xì)節(jié)信息。對(duì)低頻序列則直接采用式(1)進(jìn)行擬合，并對(duì)擬合殘差采用在線支持向量回歸(Online-SVR)建立實(shí)時(shí)殘差補(bǔ)償模型。另外，ZTD空間變化特性體現(xiàn)在均值A(chǔ)0、振幅(A1、A2)、相位(φ1、φ2)與其空間位置相關(guān)，因此建立空間位置向上述5個(gè)參數(shù)映射的SVR空間模型。本文僅以ZTD建模為例進(jìn)行說明(ZHD建模方法同ZTD)。組合模型結(jié)構(gòu)如圖2所示，各部分在下文中詳細(xì)說明。

圖2 ZTD組合預(yù)報(bào)模型框架Fig.2 Framework of ZTD combined prediction model

1.1 ZTD序列離散小波變換

離散小波變換是用一對(duì)共軛的低通和高通濾波器將信號(hào)分解成近似值(低頻)和細(xì)節(jié)值(高頻)兩部分。小波重構(gòu)為分解的逆過程，可用下式表示：

(2)

式中，a0為原始信號(hào)，ai、di分別為低頻和高頻信號(hào)。

通過實(shí)驗(yàn)發(fā)現(xiàn)， ZTD序列經(jīng)db4濾波器6層分解后的低頻序列可反映ZTD序列的主要趨勢(shì)。圖3為小波分解后的AZU1站低頻序列和高頻序列，從圖中可觀察到低頻序列曲線已變得光滑。

圖3 AZU1站ZTD采用db4小波基6層分解后的低頻序列和高頻序列Fig.3 Low-frequency and high-frequency ZTDsequences of AZU1 station generated by 6-layer decomposition with db4 wavelet basis function

1.2 高頻ZTD序列AR模型

高頻序列包含ZTD細(xì)節(jié)信息，建立AR模型的目的在于過濾白噪聲，最大程度提取出ZTD細(xì)節(jié)信息，AR模型表達(dá)式如下：

xt=φ1xt-1+φ2xt-2+…+φnxt-n+

(3)

式中，φi(i=1,2,…,n)為自回歸參數(shù)，at為白噪聲序列。

本文選用相關(guān)系數(shù)函數(shù)(ACF)和偏自相關(guān)系數(shù)函數(shù)(PACF)中截尾和拖尾進(jìn)行模型識(shí)別和定階，具體過程可參考文獻(xiàn)[7]，經(jīng)計(jì)算可對(duì)AUZ1站高頻ZTD序列建立3階AR模型(ACF拖尾，PACF截尾，模型參數(shù)為：φ1=1.145 7，φ2=-0.560 8，φ3=0.194 1)。按式(3)進(jìn)行預(yù)測(cè)，實(shí)測(cè)值與預(yù)測(cè)值對(duì)比結(jié)果見圖4。選用平均偏差(bias)和均方根誤差(RMSE)來檢驗(yàn)AR模型預(yù)測(cè)高頻ZTD序列的效果，計(jì)算公式如下：

(4)

經(jīng)計(jì)算可得bias=0.015 mm，RMSE=0.197 mm，表明預(yù)測(cè)值與實(shí)測(cè)值符合較好。結(jié)合圖4可以發(fā)現(xiàn)，預(yù)測(cè)值序列仍保留高頻振蕩特性，說明AR模型可有效提取高頻ZTD序列的細(xì)節(jié)信息。

圖4 AZU1站高頻ZTD AR模型預(yù)測(cè)值與原始值對(duì)比Fig.4 Comparison between predicted value obtained by AR model and original value of high frequency ZTD in AZU1 station

1.3 空間SVR模型

非線性規(guī)劃的支持向量回歸(SVR)問題可以表述為：給定一個(gè)訓(xùn)練樣本集T={(x1,y1),(x2,y2),…,(xl,yl)}∈(Rn×R)l，其中x∈Rn為輸入指標(biāo)向量，其分量稱為特征，y∈R為輸出指標(biāo)，在特征空間中構(gòu)造回歸函數(shù)：

f(x)=WTΦ(x)+b

(5)

式中，W為權(quán)向量，Φ(x)為核函數(shù)，可將輸入向量從Rn映射到Hilbert特征空間，b為偏置，求解帶約束的凸二次優(yōu)化問題可得W、b。

由于ZTD的空間變化特性在于均值、振幅、相位參數(shù)與其空間位置相關(guān)，統(tǒng)一以緯度B、經(jīng)度L、大地高H為SVR的輸入向量，分別對(duì)A0、A1、A2、φ1、φ2參數(shù)進(jìn)行空間回歸模型訓(xùn)練。

由SVR算法可知，核函數(shù)類型選取以及核參數(shù)γ、不敏感損失函數(shù)系數(shù)ε、懲罰參數(shù)C的取值會(huì)直接影響SVR的擬合、泛化能力。在進(jìn)行樣本訓(xùn)練時(shí)，需要進(jìn)行參數(shù)優(yōu)化。本文選取高斯徑向基核函數(shù)，參數(shù)尋優(yōu)方法為PSO(粒子群算法)，以k-fold交叉驗(yàn)證誤差作為SVR參數(shù)選擇的目標(biāo)值，詳細(xì)算法流程見文獻(xiàn)[8]。

1.4 殘差補(bǔ)償模型

Online-SVR的基本思想是通過更新回歸函數(shù)，使加入樣本后的新訓(xùn)練集繼續(xù)滿足KKT條件。Online-SVR比SVR樣本更新后的訓(xùn)練速度更快，詳細(xì)算法參照文獻(xiàn)[9-10]。

對(duì)應(yīng)一組殘差序列{v1,v2,…,vn}，Online-SVR殘差補(bǔ)償模型采用的訓(xùn)練樣本組成為：

(6)

式中，左邊矩陣為輸入向量，右邊矩陣為輸出向量。即建立滑動(dòng)時(shí)間窗口{vt-1,vt-2,…,vt-m}與vt之間的映射關(guān)系f:Rm→R。將一維時(shí)間序列轉(zhuǎn)化成矩陣形式獲得數(shù)據(jù)間的關(guān)聯(lián)關(guān)系，以挖掘到盡可能多的信息量，也稱為相空間重構(gòu)。式中m為嵌入維數(shù)，可反映轉(zhuǎn)換后矩陣蘊(yùn)涵的知識(shí)量[11]。

本文采用貝葉斯信息準(zhǔn)則(BIC)來確定維數(shù)m。由于樣本不斷更新，為提高計(jì)算效率，采用文獻(xiàn)[12]中的方法快速確定C、γ、ε(相關(guān)公式及參數(shù)詳見參考文獻(xiàn))。為提高Online-SVR殘差補(bǔ)償模型的訓(xùn)練速度，保證訓(xùn)練樣本之間的相關(guān)性，樣本容量不宜過大。本文初步確定樣本容量閾值為10，當(dāng)達(dá)到閾值后，每加入1個(gè)新樣本，就需要移除最早時(shí)刻的1個(gè)樣本。Online-SVR殘差補(bǔ)償模型預(yù)報(bào)流程如圖5所示。

圖5 Online-SVR殘差補(bǔ)償模型預(yù)報(bào)流程Fig.5 Prediction process of online-SVRresidual compensation model

2 模型精度檢驗(yàn)

采用2016～2017年ZTD和ZHD序列數(shù)據(jù)進(jìn)行建模，以2018年ZTD和ZHD序列為真值進(jìn)行預(yù)報(bào)精度檢驗(yàn)，模型精度評(píng)價(jià)指標(biāo)選用bias和RMSE。由于建模采用多個(gè)GNSS基站數(shù)據(jù)，ZTD/ZHD預(yù)報(bào)模型的AR階數(shù)和殘差補(bǔ)償模型的嵌入維數(shù)采用建模站點(diǎn)的平均值序列進(jìn)行計(jì)算。表1為最后確定的ZTD/ZHD模型參數(shù)。

表1 模型參數(shù)設(shè)置

根據(jù)建立的區(qū)域?qū)α鲗友舆t預(yù)報(bào)模型可計(jì)算區(qū)域內(nèi)任意位置對(duì)應(yīng)某個(gè)年積日的ZTD/ZHD，計(jì)算流程如下：

1) 輸入B、L、H，由空間SVR模型計(jì)算A0、A1、A2、φ1、φ2；

2) 輸入doy，按式(1)計(jì)算ZTD/ZHD初始值；

3) 殘差補(bǔ)償模型補(bǔ)償ZTD/ZHD初始值；

4) 按式(2)疊加AR模型預(yù)報(bào)的高頻ZTD/ZHD值；

5) 輸出ZTD/ZHD最終預(yù)報(bào)值。

任意位置步驟3)、4)數(shù)值可由周邊一定數(shù)量建模站點(diǎn)的預(yù)報(bào)值反距離加權(quán)內(nèi)插得到。

2.1 內(nèi)符合精度檢驗(yàn)

從建模站點(diǎn)中挑選觀測(cè)數(shù)據(jù)較完整的80個(gè)基站2018年的ZTD/ZHD值與模型預(yù)報(bào)值進(jìn)行對(duì)比，檢驗(yàn)?zāi)Ｐ偷目煽砍潭取?nèi)符合精度統(tǒng)計(jì)如圖6所示。

圖6 建模站點(diǎn)ZTD/ZHD內(nèi)符合檢驗(yàn)精度統(tǒng)計(jì)Fig.6 Statistics of internal coincidence accuracy of ZTD/ZHD of modeling stations

統(tǒng)計(jì)數(shù)據(jù)表明，ZHD和ZTD的平均bias均為0.3 mm左右，RMSE分別為5.67 mm和2.34 cm。模型預(yù)報(bào)ZTD精度較ZHD低，原因在于ZTD包含變化相對(duì)更劇烈的ZWD。但總體上模型內(nèi)符合程度較好。

2.2 外符合精度檢驗(yàn)

將24個(gè)測(cè)站2018年的ZTD/ZHD值與模型預(yù)測(cè)值進(jìn)行對(duì)比，測(cè)試模型的適用性(泛化能力)。圖7為模型外符合精度，從圖中可以看出，外符合精度稍低于內(nèi)符合精度(mm級(jí)差別)，說明模型泛化能力好。此外，ZTD和ZHD的RMSE最大、最小值相差較小，表明模型在區(qū)域內(nèi)的預(yù)報(bào)精度較均勻。

3 模型應(yīng)用效果測(cè)試

本文將該模型應(yīng)用到偽距單點(diǎn)定位中，通過分析模型對(duì)定位結(jié)果產(chǎn)生的影響來檢驗(yàn)?zāi)Ｐ皖A(yù)測(cè)ZTD/ZHD的可靠性。實(shí)驗(yàn)數(shù)據(jù)為測(cè)站(BRIB、FARB、MASW、MLFP、OHLN)2018-01-22全天的原始觀測(cè)值。對(duì)流層延遲分別采用模型預(yù)報(bào)值改正和不加改正兩種模式。

3.1 偽距單點(diǎn)定位解算策略

觀測(cè)值為無電離層組合偽距(削弱電離層延遲)，截止高度角取15°，采用精密星歷，對(duì)流層延遲映射函數(shù)為VMF1。為盡可能削弱其他誤差對(duì)定位結(jié)果的影響，本文實(shí)驗(yàn)還進(jìn)行地球自轉(zhuǎn)改正、衛(wèi)星天線相位中心改正、相位纏繞改正和DCB補(bǔ)償。

3.2 偽距單點(diǎn)定位結(jié)果分析

兩種模式下偽距單點(diǎn)定位在N、E、U方向的偏差如表2所示。經(jīng)模型ZTD預(yù)測(cè)值改正后U方向的定位偏差從6～8 m縮小至dm級(jí)，同文獻(xiàn)[13]中對(duì)流層延遲對(duì)單點(diǎn)定位U方向影響最大(可達(dá)7～15 m)的結(jié)論一致。另外，本實(shí)驗(yàn)取截止高度角15°使得低高度角衛(wèi)星不參與解算，對(duì)流層延遲對(duì)平面定位結(jié)果的影響會(huì)進(jìn)一步降低。因此改正后N、E方向的偏差與不改正模式差別較小。

表2 兩種模式下偽距單點(diǎn)定位N、E、U方向偏差對(duì)比

表3為兩種模式下偽距單點(diǎn)定位在N、E、U三個(gè)方向的中誤差。從表中可以看出，二者基本保持一致，說明本文所建的組合模型在對(duì)偽距觀測(cè)值的對(duì)流層延遲進(jìn)行有效改正時(shí)，未引入多余噪聲而降低單點(diǎn)定位精度。

表3 兩種模式下偽距單點(diǎn)定位N、E、U方向中誤差對(duì)比

4 結(jié) 語

本文提出一種結(jié)合小波變換、傅里葉級(jí)數(shù)擬合、AR和SVR的對(duì)流層延遲建模新方法，選用NOAA CORS網(wǎng)在加州境內(nèi)的94個(gè)GNSS基站2016～2017年數(shù)據(jù)建立區(qū)域ZTD/ZHD數(shù)值預(yù)報(bào)模型，并對(duì)模型內(nèi)、外符合精度進(jìn)行檢驗(yàn)。總體來說，預(yù)報(bào)模型內(nèi)、外符合精度相差較小，表現(xiàn)出較好的泛化性能。

分析組合模型結(jié)構(gòu)可以看出，該模型包含函數(shù)和預(yù)測(cè)兩部分，使得模型預(yù)測(cè)值既符合ZTD/ZHD變化趨勢(shì)又可保留部分細(xì)節(jié)變化屬性，能更好地描述ZTD/ZHD的時(shí)空變化特征。24個(gè)測(cè)站1 a的ZTD/ZHD預(yù)測(cè)結(jié)果表明，平均bias為-2.02 mm/-0.98 mm，RMSE為3.07 cm/6.10 mm，精度優(yōu)于現(xiàn)有的大部分區(qū)域?qū)α鲗友舆t模型。將該模型應(yīng)用到偽距單點(diǎn)定位中可取得較好的對(duì)流層延遲改正效果，可顯著提高單點(diǎn)定位U方向精度。

綜上所述，該組合模型具有較高的預(yù)報(bào)精度和可靠性，在GNSS領(lǐng)域具有一定的應(yīng)用價(jià)值。

致謝：感謝NOAA CORS 提供的部分GNSS基站觀測(cè)數(shù)據(jù)。