GNSS觀測數據ZTD建模的質量控制方法

丁君生，陳俊平，王君剛

(1. 中國科學院上海天文臺，上海 200030； 2. 中國科學院大學天文與空間科學學院，北京 100049； 3. 德國地學中心，波茲坦14473)

0 引 言

電磁波信號在穿過大氣層時受到氣體分子的影響，被延遲和彎曲，從而使得測量距離產生偏差[1]。其中穿過電離層產生的偏差叫做電離層延遲，而穿過未被電離的中性大氣(包含對流層和平流層)產生的偏差統稱為對流層延遲[2]。電離層延遲可以通過雙頻差分技術基本消除，而對流層具有非色散性質，其延遲無法通過雙頻數據改正[3]。在全球衛星導航系統(Global navigation satellite system，GNSS)數據處理中，對流層延遲通常作為待估參數解算或通過模型進行改正[4]。

模型改正法可以分為兩類，第一類是基于實測氣象參數的對流層模型，如Saastamoinen[5](以下簡稱為Saas)和Hopfield[6]等，其改正精度可達厘米級或分米級[7]。由于這類模型高度依賴實測氣象參數，且有學者發現這類模型計算的ZTD與經驗模型相比精度上沒有優勢[8]，故而這類模型在實時導航中的應用偏少。第二類是基于對流層關鍵參量的經驗模型，這類模型不依賴實測氣象參數，使用方便，僅需提供時間以及測站信息，便能得到厘米級精度的ZTD。其中有基于標準大氣模型參數的UNB系列模型[9]和EGNOS模型[10]，UNB3模型在北美區域估計的ZTD平均誤差約為2 cm[11]，EGNOS總體精度與采用實測氣象參數的Saas和Hopfield模型相當[8]。還有基于再分析資料的GPT系列[12-14]、TropGrid系列[15-16]和IGGtrop系列模型[17]，GPT模型基于歐洲中期天氣預報中心(European Centre for Medium-Range Weather Forecasts，ECMWF)提供的數值天氣模型(Numerical Weather Model，NWM)產品，可提供測站氣象數據，然后采用Saas等模型計算ZTD，精度優于EGNOS、UNB3模型。TropGrid基于美國國家環境預報中心(National Centers for Environment Prediction，NCEP)的NWM產品建立，與EGNOS相比全球平均精度提高了25%[18]。IGGtrop模型基于美國國家環境預報中心(National Centers for Environment Prediction，NCEP)再分析資料建立，精度優于EGNOS、UNB3、UNB3m[17]。

參數估計法將ZTD作為未知參數參與GNSS定位解算，其將ZTD的主項，即天頂對流層靜力學延遲(Zenith hydrostatic delay，ZHD)作為延遲的先驗值，將天頂對流層濕延遲(Zenith wet delay，ZWD)作為ZHD的改正值[19]。參數估計法精度最高，可達毫米級[20]，但是受限于GNSS測站的分布，無法獲得任意位置的對流層延遲估計。

隨著GNSS站網的加密以及數據傳輸速率的提升，高空間分辨率的ZTD數據出現爆炸式地增長，部分學者開始采用一些機構解算的GNSS對流層延遲數據，建立對流層延遲模型。這類模型保留了經驗模型不依賴實測氣象參數、建模簡單、使用方便和精度優良的優點，同時能夠最大程度保證在GNSS應用中的自洽。文獻[21]采用國際GNSS服務組織(International GNSS Survice，IGS)提供的對流層天頂延遲產品，建立了一種簡單的全球對流層天頂延遲模型，全球范圍內RMS為4.9 cm。受限于IGS測站數量不足，該模型并不精細。文獻[4]利用上海天文臺GNSS分析中心解算的陸態網對流層數據建立了中國區域對流層天頂延遲網格模型SHAtrop，該模型建模精細，精度良好。然而在整個模型的建立過程中，缺乏有效的質量控制手段，不能夠給其他區域建立類似的模型提供很好的參考。同時，該模型在ZTD時間序列擬合中，采用了帶相位的周期函數，使用該函數擬合出的初相并不穩定，可能在制作網格以及使用網格的插值過程中引入誤差。

為了填補這些內容的空缺，本文提出了一套ZTD建模質量控制方法，采用內華達大地測量實驗室(Nevada Geodetic Laboratory，NGL)解算的高空間分辨率GNSS對流層數據，選取了近10年德國及周邊區域[47°N-55°N，5°E-15°E]183個測站的實測ZTD，對該方法進行了檢驗。該方法從建模數據量的選擇、網格分辨率對模型精度的影響和模型產品穩定度等方面對模型進行了質量控制。同時，本文提出了無相位周期擬合模型并對比分析了其相對于文獻[4]中含相位模型的優勢。最后，本文建立了在該質量控制方法和新擬合方法下所選區域的ZTD模型。新模型穩定可靠，精度優良，檢驗了該質量控制方法的可行性。

1 數據處理策略與模型質量控制

1.1 建模步驟

相比于傳統模型，基于GNSS觀測數據的ZTD模型與氣象參數無關，其建模數據來源于區域或全球分布的GNSS永久觀測站的實測ZTD，模型建立過程中充分考慮ZTD的時空分布特性，能夠最大程度地保證模型在GNSS應用中的自洽。文中，為了使得模型網格盡可能的平滑以降低插值帶來的誤差，在數據處理策略上，首先剝離高程給ZTD帶來的影響，即將測站ZTD歸化到橢球面上，然后再網格化。

本文的建模步驟為：1)獲取區域的GNSS測站坐標信息以及長時間的ZTD數據；2)分析測站高程與ZTD年均值的關系，采用指數模型計算ZTD隨高程的衰減系數，根據計算獲取的系數將ZTD歸化到橢球面上；3)采用頻譜分析方法提取測站ZTD時間序列里的周期信號，使用周期函數對時間序列進行擬合，獲取各個測站的擬合參數；4)分析各個參數水平方向的分布特點，對各個測站的擬合參數采用雙線性插值方法在經緯度上進行網格化，獲取各網格點的參數，制作成網格文件。

模型的使用過程與建立過程互為逆過程，其步驟為：1)獲取網格文件，通過搜索和用戶最近的4個網格點，使用雙線性插值方法獲取用戶所在經緯度的周期函數擬合參數；2)將參數帶入周期擬合模型，獲取給定年積日下用戶所在經緯度橢球面上的ZTD；3)使用指數模型，帶入衰減系數參數，將橢球面上的ZTD改正到用戶所在高程，求得測站的ZTD。

建模及使用具體流程如圖1所示，圖1上部點劃線框中為模型建模過程，下部點劃線框中為模型的使用過程。對于用戶而言，僅需提供測站坐標和年積日，代入模型，便能獲取實時ZTD。

1.2 數據質量評估

高空間分辨率以及高精度的實測GNSS ZTD是模型建立的基礎。2017年11月5日，NGL的Geoff Blewitt等人開放提供了自1996年以來全球超過1.6萬個站點超過3400萬天次的對流層產品[23]，每年新增測站數量約為1000個。NGL提供的對流層產品使用了JPL/Caltech提供的軌道和鐘差，使用GIPSY軟件進行解算，時間分辨率為5 min，產品遵循IGS標準格式。相比于IGS對流層產品，NGL對流層產品測站分布范圍更廣，分布密度更大，可用時間更長。NGL對流層數據的開放對于對流層延遲的精細建模研究具有重要的意義。

相比于全球模型，區域模型由于空間跨度較小，建模更加精細，往往能夠獲得更高的精度。鑒于德國區域NGL測站數量多且分布均勻，且該區域是用于建立歐洲氣象基礎模型的主要觀測區域之一，為便于與GPT2w等基于NWM再分析產品的模型進行對比分析，本文選取了2009年1月1日至2018年12月31日期間德國及周邊區域[47°N-55°N，5°E-15°E]217個測站的ZTD數據，用于實驗分析。

NGL并未公布其對流層產品的精度，為了評估本文使用的NGL對流層產品的精度，以下以IGS ZTD為參考，統計了10年期間所選測區的10個共有測站的NGL ZTD的平均偏差(BIAS)和均方根誤差(RMS)。統計方法見式(1)，統計結果見表1，從表1可以看出，測區內的NGL ZTD和IGS ZTD一致性良好，10個測站BIAS均值為0.87 mm，RMS均值為3.88 mm。該結果表明NGL解算的ZTD具有和IGS解算的ZTD相同的精度，可用于對流層延遲精細建模。(IGS提供的最終ZTD產品精度為4 mm；http://www.igs.org/products，2019.7)

(1)

式中：B為BIAS，R為RMS，ZIGS為IGS ZTD序列，ZNGL為NGL ZTD序列，N為序列長度。

1.3 數據量對模型的影響

圖2中柱狀圖為所選測區內各個測站數據可用時間的頻數分布直方圖，從圖中可知各個測站之間可用數據差異較大。前人的研究中[4,21]，僅簡單剔除了數據量偏少的測站，未在數據量的選擇對建模精度的影響上進行深入研究。為了解決不同數據量的測站用來一起建模是否合理、數據量越多是否越有益于提高模型精度等質量控制問題，實驗選取了數據量滿10年的15個測站，計算了分別使用1年至10年數據量建模，測站的平均擬合RMS值，結果如圖2所示?？梢园l現，隨著數據量的增多，RMS值在增大，幅度約為0.8 mm/年，但在數據量超過4年后，RMS逐漸平緩。由此可以得到結論：1)ZTD在年與年之間存在差異，這種差異在數據量超過1年的情況下會使得RMS值變大；2)這種差異的周期約為4年，在數據量超過4年后，RMS值趨于平緩，因此在數據量的選取上，建議不超過4年；3)這種差異的幅度為毫米級，將不同數據量的測站一起建模對精度影響有限。在本文中，考慮到在時間序列擬合過程中使用了年周期擬合模型，實驗對217個測站進行篩選，在剔除數據量不足1年的測站后，選取了2015-2018共4年期間183個測站的ZTD數據進行建模。

圖2 數據可用時間頻數分布和建模數據量 與模型RMS之間的關系Fig.2 Data available time frequency distribution & Relationship between modeling data volume and model RMS

表1 NGL解算的德國區域IGS測站ZTD精度Table 1 ZTD accuracy of the IGS station in The German region calculated by NGL (mm)

1.4 網格分辨率對精度的影響

測站空間分辨率的大小直接決定了模型的精細程度，而網格分辨率的選擇與測站空間分辨率相關。前人的研究[4]中未闡明其網格分辨率的選擇依據，無法獲知網格分辨率對模型精度的影響。為了探究網格分辨率對模型的影響，以確定合適的網格分辨率，本文統計了在使用相同的數據、相同的建模方式下劃分成不同大小網格下的模型精度(BIAS和RMS)以及模型使用的成本(網格文件大小及計算耗時)。實驗結果如表2所示，由于網格數量與網格點分布相關，難以統計，本文中以網格點數量代替，表中平均測站數為測區內用以建模的測站個數與網格點個數之商。從表2可以發現：1)隨著網格密度的增大，網格點數量、網格文件大小呈指數增長，計算耗時前期增長不大，后期增長迅速；2)RMS隨著網格的變小逐漸變小，當分辨率小于1°×1°后，RMS趨于穩定，BIAS隨著網格的變小有明顯的下降；3)網格分辨率在小于1°×1°之后，網格平均測站數小于1。綜合考慮，1°×1°的網格分辨率在本文中為最佳網格劃分標準，在該標準下，既保證了高精度(BIAS小于1 mm)和高效率(耗時較低)，又能保證了低成本(網格文件占用內存較小)。從以上結果可以得到網格分辨率選取的一般標準：在保證平均每個網格內至少1個測站的前提下選擇大網格。在該標準下，能夠在充分發揮測站分辨率大的優勢同時控制建模及使用成本。

1.5 模型穩定性

由1.3節可知，建模數據時長的增加并不能有效地提高模型精度，可見ZTD具有較高的穩定性。為了檢驗該結論并探究模型精度在模型使用過程中隨時間的變化，實驗選取了數據量滿10年的15個測站2009-2012年共4年的數據進行建模，然后對2013-2018共6年時間的ZTD進行預報。圖3為15個測站的預報殘差平均值的時間序列，圖中點劃線為殘差序列的線性擬合，可以發現，在長時間的預報過程中殘差序列存在逐漸變大的趨勢，殘差序列變化率為1.42 mm/年?？紤]到模型精度為厘米級，該趨勢變化并不大，模型精度整體相對穩定，無需頻繁(每年)更新，該模型可以在建模后的較長一段時間以其標稱精度提供服務。

圖3 模型預報殘差Fig.3 Model prediction residual

2 結果與討論

為了檢驗上述質量控制方法下的模型精度，本文對在該區域選取的183個測站數據，進行了ZTD時空特性分析，建立了相應的ZTD模型，該模型以下稱為SHAtropDE(DE為德國的簡稱)。

2.1 高程歸化

圖4為各測站高程以及各測站的ZTD均值，本文中的測站高程均為基于橢球面的大地高。從圖中可以發現，德國地勢整體較為平緩，南部由于靠近阿爾卑斯山脈高程出現陡增，而ZTD年均值正好相反，在南部出現陡降，測站高程與ZTD近似呈反比關系，即測站ZTD年均值隨著測站高程增加減小。ZTD與測站高程的關系通常用指數模型或者線性模型來表示[22]，以下為指數模型關系式：

Z(h)=Z0×exp(βh)

(2)

式中：Z(h)為測站所在高程ZTD，Z0為測站在橢球面上的ZTD，h為測站高程(m)，β為ZTD隨高程的衰減系數。

實驗通過對測站4年的ZTD年均值與測站高程進行指數模型擬合，得衰減系數β=-1.24×10-4。獲取衰減系數后，即可將ZTD歸化到橢球面上。

表2 不同網格分辨率下模型精度及使用成本Table 2 Model accuracy and cost of use at different grid resolutions

圖4 測站高程和測站ZTD年均值Fig.4 Station elevation and station ZTD annual mean

2.2 時間序列周期擬合

由2.1節可知測站ZTD與測站高程相關，而對于單個測站而言，其高程固定，ZTD隨時間的變化規律與高程無關。圖5為PTBB測站的時間序列，可以發現ZTD時間序列存在明顯的周期信號。時間序列的周期信號通常通過傅里葉頻譜分析[24]來獲取，實驗對全部測站數據做傅里葉頻譜分析后發現全部測站均存在顯著的年周期信號。

在處理ZTD隨時間的變化中，文獻[4]采用了年周期+半年周期擬合模型，公式見式(3)。該模型含有相位參數，為了克服由該相位參數帶來的模型誤差，本文提出使用無相位擬合模型。為了判斷在本實驗所在區域中，年周期+半年周期是否同樣優于僅年周期模型，本文提出了兩種無相位模型，第一種為無相位僅年周期模型，見式(4)；第二種為無相位年周期+半年周期模型，見式(5)。

(3)

(4)

(5)

式中：A,A1,B1為年周期項振幅，B,A2,B2為半年周期項振幅，P1,P2分別為年周期和半年周期項相位，C為常數項，τ為年積日。

實驗對式(4)和式(5)模型、式(3)和式(5)模型下的RMS進行了對比實驗。實驗結果如圖6所示。圖中分別以“H”“N”“O”代表無相位僅年周期模型、無相位年周期+半年周期模型和含相位年周期+半年周期模型。

從圖6 (a)子圖可以發現，無相位僅年周期模型和無相位年周期+半年周期模型RMS序列接近，但RMS差值序列全部為正，說明在本實驗所在區域內所有測站在附加半年周期模型下精度高于僅年周期模型。根據圖6 (b)子圖，使用無相位模型精度高于含相位模型，證明了新模型的優越性。PTBB站使用式(5)模型的擬合結果以及擬合殘差也繪制在圖5中，圖中對殘差取了絕對值，從圖5可以發現，擬合殘差基本控制在±10 cm以內，且殘差序列呈現噪聲特性，不存在殘余的周期信號，擬合情況良好。

圖5 測站ZTD時間序列與周期函數擬合結果及擬合殘差Fig.5 Station ZTD time series and cycle function fit results and fitting residuals

圖6 三種擬合模型下各個測站RMS比較Fig.6 RMS comparison of each station under three models

2.3 水平方向網格化

實驗選取183個測站中的170個用于建模，剩余的13個用于檢驗。圖7為建模站做年周期+半年周期擬合的5個擬合函數參數和擬合殘差的中誤差。從圖中可以發現：年周期項振幅A1表現為西北部大、東南部小，并有向東南方向層進減小的趨勢，而B1無顯著地域特性；半年項振幅A2可以看出存在南大北小的地域分布，而B2同樣無明顯地域特征；常數項C在西南區域較大，在南邊界靠近阿爾卑斯山處較小，其余區域呈現由北向南增大的特點；擬合殘差中誤差RMS呈現從臨海的西北區域向靠近阿爾卑斯山脈的東南內陸逐漸減小的特點，可見水汽活躍程度極大影響擬合精度。測站中擬合殘差中誤差最大為3.71 cm(VLIE∶5.09E, 53.30 N)，最小為2.90 cm(PAT2∶11.46E, 47.21 N)。

2.4 模型精度評估

為了評估SHAtropDE的精度，本文利用SHAtropDE計算了所選測區內183個測站2015-2018共4年的ZTD，同時計算了EGNOS、UNB3m和GPT2w模型在相同測站、相同時間段的ZTD，并以NGL解算的ZTD作為參考值，計算了各個模型的BIAS和RMS，計算公式見式(1)。各模型下各個測站的BIAS如圖8(a)至(d)，RMS如圖8(e)至(h)。

從BIAS和RMS可以看出，UNB3m模型整體精度最低，EGNOS和GPT2w精度相當，SHAtropDE模型精度最高。EGNOS和GPT2w的BIAS以及SHAtropDE的RMS均呈現南高北低的地域特性，SHAtropDE模型BIAS最小且分布均勻。

表3統計了170個建模站和13個檢驗站在四個模型下的BIAS和RMS信息(平均值，最小值，最大值)，從表3可以得到以下結論：SHAtropDE模型在建模站和驗證站精度相當，平均BIAS小于1 mm，平均RMS為3.4 cm；SHAtropDE模型的整體精度最高，建模站相對于UNB3m、EGNOS、GPT2w+Saas平均改善了42.4%、35.8%、33.3%；驗證站相對于UNB3m、EGNOS、GPT2w+Saas平均改善了30.8%、40.7%、38.6%。SHAtropDE模型基于實測ZTD建立，相對于傳統經驗模型有良好的精度改善，能夠滿足德國區域GNSS用戶高精度ZTD改正需求。

3 結 論

本文提出了一套ZTD建模質量控制方法，并采用NGL解算的高空間分辨率GNSS對流層數據，選取了近10年德國及周邊區域[47°N-55°N，5°E-15°E]183個測站的實測ZTD，對該方法進行了驗證，得出了以下結論：

1)ZTD在年與年之間存在大小約為0.8 mm/年的擬合差異，這種差異的周期約為4年，使用超過4年的數據對模型精度提升作用不明顯；

2)網格分辨率的選取應以保證平均每個網格內至少一個測站為標準。在該標準下，能夠充分發揮測站分辨率大優勢并控制建模及使用成本；

圖7 各測站擬合殘差及擬合參數Fig.7 Fitting residuals and fitting parameters of each station

圖8 各測站不同ZTD模型的系統偏差(cm)和RMS精度統計(cm)Fig.8 System deviation (cm) and RMS accuracy statistics (cm) for different ZTD models at each station

表3 SHAtropDE模型的精度與其他模型的對比Table 3 Comparison of the accuracy of the SHAtropDE model with other models (cm)

3)模型在長時間的使用過程中存在精度略微退化的現象，殘差序列變化率為1.42 mm/年。模型精度整體相對穩定，無需頻繁(每年)更新，可以在建模后的較長一段時間以其標稱精度提供服務。

在本文質量控制方法下建立的新模型平均RMS為3.4 cm，相對于UNB3m、EGNOS、GPT2w+Saas平均改善了42.4%、35.8%、33.3%，能夠滿足德國區域GNSS用戶高精度實時ZTD改正需求。該方法為此類ZTD模型的建立制定了技術規范，能夠給其他區域的建立類似模型提供良好的參考作用，對于ZTD建模研究具有一定的參考價值。