一種區域精密天頂對流層延遲組合預報模型

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

1 甘肅農業職業技術學院，蘭州市段家灘路425號，730030 2 中國地震局地質研究所地震動力學國家重點實驗室，北京市華嚴里甲1號，100029

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

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

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

1 建模方法及原理

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

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

為研究ZTD的時間變化特征，隨機選取1個站點(AZU1)2016～2018年的ZTD時間序列進行快速傅里葉變換，圖1為時間序列和對應的頻譜圖。

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

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

(1)

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

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

圖2 ZTD組合預報模型框架Fig.2 Framework of ZTD combined prediction model

1.1 ZTD序列離散小波變換

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

(2)

式中，a0為原始信號，ai、di分別為低頻和高頻信號。

通過實驗發現， ZTD序列經db4濾波器6層分解后的低頻序列可反映ZTD序列的主要趨勢。圖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細節信息，建立AR模型的目的在于過濾白噪聲，最大程度提取出ZTD細節信息，AR模型表達式如下：

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

(3)

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

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

(4)

經計算可得bias=0.015 mm，RMSE=0.197 mm，表明預測值與實測值符合較好。結合圖4可以發現，預測值序列仍保留高頻振蕩特性，說明AR模型可有效提取高頻ZTD序列的細節信息。

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

1.3 空間SVR模型

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

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

(5)

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

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

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

1.4 殘差補償模型

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

對應一組殘差序列{v1,v2,…,vn}，Online-SVR殘差補償模型采用的訓練樣本組成為：

(6)

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

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

圖5 Online-SVR殘差補償模型預報流程Fig.5 Prediction process of online-SVRresidual compensation model

2 模型精度檢驗

采用2016～2017年ZTD和ZHD序列數據進行建模，以2018年ZTD和ZHD序列為真值進行預報精度檢驗，模型精度評價指標選用bias和RMSE。由于建模采用多個GNSS基站數據，ZTD/ZHD預報模型的AR階數和殘差補償模型的嵌入維數采用建模站點的平均值序列進行計算。表1為最后確定的ZTD/ZHD模型參數。

表1 模型參數設置

根據建立的區域對流層延遲預報模型可計算區域內任意位置對應某個年積日的ZTD/ZHD，計算流程如下：

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

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

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

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

5) 輸出ZTD/ZHD最終預報值。

任意位置步驟3)、4)數值可由周邊一定數量建模站點的預報值反距離加權內插得到。

2.1 內符合精度檢驗

從建模站點中挑選觀測數據較完整的80個基站2018年的ZTD/ZHD值與模型預報值進行對比，檢驗模型的可靠程度。內符合精度統計如圖6所示。

圖6 建模站點ZTD/ZHD內符合檢驗精度統計Fig.6 Statistics of internal coincidence accuracy of ZTD/ZHD of modeling stations

統計數據表明，ZHD和ZTD的平均bias均為0.3 mm左右，RMSE分別為5.67 mm和2.34 cm。模型預報ZTD精度較ZHD低，原因在于ZTD包含變化相對更劇烈的ZWD。但總體上模型內符合程度較好。

2.2 外符合精度檢驗

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

3 模型應用效果測試

本文將該模型應用到偽距單點定位中，通過分析模型對定位結果產生的影響來檢驗模型預測ZTD/ZHD的可靠性。實驗數據為測站(BRIB、FARB、MASW、MLFP、OHLN)2018-01-22全天的原始觀測值。對流層延遲分別采用模型預報值改正和不加改正兩種模式。

3.1 偽距單點定位解算策略

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

3.2 偽距單點定位結果分析

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

表2 兩種模式下偽距單點定位N、E、U方向偏差對比

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

表3 兩種模式下偽距單點定位N、E、U方向中誤差對比

4 結 語

本文提出一種結合小波變換、傅里葉級數擬合、AR和SVR的對流層延遲建模新方法，選用NOAA CORS網在加州境內的94個GNSS基站2016～2017年數據建立區域ZTD/ZHD數值預報模型，并對模型內、外符合精度進行檢驗。總體來說，預報模型內、外符合精度相差較小，表現出較好的泛化性能。

分析組合模型結構可以看出，該模型包含函數和預測兩部分，使得模型預測值既符合ZTD/ZHD變化趨勢又可保留部分細節變化屬性，能更好地描述ZTD/ZHD的時空變化特征。24個測站1 a的ZTD/ZHD預測結果表明，平均bias為-2.02 mm/-0.98 mm，RMSE為3.07 cm/6.10 mm，精度優于現有的大部分區域對流層延遲模型。將該模型應用到偽距單點定位中可取得較好的對流層延遲改正效果，可顯著提高單點定位U方向精度。

綜上所述，該組合模型具有較高的預報精度和可靠性，在GNSS領域具有一定的應用價值。

致謝：感謝NOAA CORS 提供的部分GNSS基站觀測數據。