張繼洪,謝劭峰,張亞博,曾 印,唐友兵,熊 思

(1.桂林理工大學(xué) 測繪地理信息學(xué)院，廣西 桂林 541004；2.湖北科技學(xué)院 資源環(huán)境科學(xué)與工程學(xué)院，湖北 咸寧 437100)

0 引言

對流層延遲是電磁波信號在傳播過程中的主要誤差源，其在天頂方向上的延遲約為2 m，隨著衛(wèi)星高度角的降低，其影響將增至20 m[1-2]。對流層天頂延遲(Zenith Total Delay，ZTD)模型對全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System，GNSS)實時高精度導(dǎo)航定位及GNSS水汽反演具有重要意義。在已有的ZTD模型中，Hopfield模型[3]和Saastamoinen模型[4]雖然可以提供精度較高的ZTD，但其應(yīng)用時需要輸入實測氣象參數(shù)；UNB3m模型[5]和EGNOS模型[6]雖然無需輸入實測氣象參數(shù)，但其較低的時空分辨率影響了ZTD計算精度。近年來，眾多學(xué)者利用大氣再分析資料計算的ZTD進行ZTD模型的建立，但由于ZTD在垂直方向的變化比水平方向的變化劇烈且復(fù)雜，大氣再分析資料格網(wǎng)點與用戶位置存在高度差，由格網(wǎng)點ZTD插值到用戶位置時將產(chǎn)生較大誤差。為此，有學(xué)者在對流層垂直剖面函數(shù)方面開展了廣泛的研究[7-8]，通常采用二次多項式[9]、負指數(shù)函數(shù)[10-15]和高斯型函數(shù)[16]進行模型構(gòu)建。研究發(fā)現(xiàn)，利用高斯型函數(shù)擬合垂直剖面ZTD效果略優(yōu)于二次多項式和負指數(shù)型函數(shù)。以上模型還存在時空分辨率較低、在部分區(qū)域(如廣西地區(qū))適應(yīng)性較差等不足。因此亟需構(gòu)建適用于廣西地區(qū)高精度、高時空分辨率的ZTD垂直剖面模型。

本文利用廣西地區(qū)2016—2017年高時空分辨率的ERA5再分析資料積分計算的分層ZTD，經(jīng)高斯函數(shù)對垂直剖面ZTD進行模型化，分析模型系數(shù)b與模型系數(shù)c的精細時空特性，為廣西地區(qū)高精度ZTD垂直剖面模型的構(gòu)建提供參考。

1 數(shù)據(jù)來源與計算原理

1.1 數(shù)據(jù)來源

ERA5大氣再分析資料(https:∥apps.ecmwf.int/datasets/data/interim-full-daily)是歐洲中期天氣預(yù)報中心(European Centre for Medium-Range Weather Forecasts,ECMWF)第5代全球氣候再分析數(shù)據(jù)集，該數(shù)據(jù)集提供逐小時(UTC時)覆蓋全球的氣壓、溫度和比濕等氣象參數(shù)，包括地表數(shù)據(jù)和分層數(shù)據(jù)，水平分辨率為0.25°×0.25°(緯度×經(jīng)度)，垂直方向包括37層，以標準氣壓層進行分層。ERA5大氣再分析資料是目前時空分辨率最高的大氣再分析資料，是進行大氣科學(xué)研究和對流層模型構(gòu)建的重要數(shù)據(jù)源。文獻[17]利用桂林市GNSS觀測站數(shù)據(jù)解算的ZTD對ERA5大氣再分析資料在桂林市計算ZTD的精度進行了評估，發(fā)現(xiàn)ERA5大氣再分析資料在該地區(qū)計算ZTD具有較高的精度。本文選取2016—2017年廣西地區(qū)(104°～113°E，20°～27°N)的ERA5大氣再分析資料。

1.2 計算原理

根據(jù)ERA5大氣再分析資料提供的地表和分層數(shù)據(jù)，利用積分法計算廣西地區(qū)每個格網(wǎng)點的垂直分層ZTD[18]：

(1)

(2)

(3)

式中，e為水汽壓，單位hPa；Sh為比濕；P為氣壓，單位hPa；N為總折射率；T為溫度，單位K；H為高程，單位km；hlow為大氣積分計算的最底層高度；htop為大氣積分計算的最頂層高度；k1，k2，k3均為常系數(shù)，其值分別為77.604 K/Pa，64.79 K/Pa，375 463 K2/Pa。

由于大氣再分析資料ERA5頂層存在殘余大氣，且主要存在對流層天頂靜力學(xué)延遲(Zenith Hydrostatic Delay，ZHD)，而對流層濕延遲(Zenith Wet Delay，ZWD)較小，可忽略不計。為了提高ZTD的積分計算精度，利用Saastamoinen模型計算ERA5再分析資料頂層上大氣殘余ZHD，將其附加在每一層的積分結(jié)果上：

(4)

式中，Ptop為對流層頂層大氣壓，單位hPa；φ為緯度；ZHDSaas為Saastamoinen模型計算的ERA5層頂殘余ZHD。由于ERA5數(shù)據(jù)是離散的分層數(shù)據(jù)而非連續(xù)數(shù)據(jù)，因此將式(3)進行離散化：

(5)

式中，Ni為第i層大氣的總折射率；ΔHi為第i層大氣的厚度；n為積分層數(shù)。

鑒于高斯函數(shù)對垂直剖面ZTD的擬合效果略優(yōu)于二次多項式和指數(shù)形式，因此，本文選用高斯函數(shù)對垂直剖面ZTD進行擬合：

(6)

式中，a為模型系數(shù)，單位m；b，c也為模型系數(shù)，單位km；e為自然常數(shù)；h為高程，單位km；ZTD(h)為高程h處的ZTD值，單位m。為了驗證高斯函數(shù)對垂直剖面ZTD的擬合效果，選取廣西地區(qū)具有代表性的4個格網(wǎng)點，通過積分法計算2017年1月1日00:00(UTC時)的分層ZTD，采用高斯型函數(shù)對其進行擬合，結(jié)果如圖1和表1所示。

(a)26°N,111°E

表1 ZTD在垂直剖面的高斯函數(shù)擬合偏差和均方根誤差

由圖1可知，每個格網(wǎng)點不同高程的ZTD散點與高斯函數(shù)擬合線均具有很高的吻合度，4個格網(wǎng)點在10 km以下均表現(xiàn)出較大的偏差，5 km以下表現(xiàn)出正偏差，5～10 km表現(xiàn)出負偏差，4個格網(wǎng)點在20 km以上均表現(xiàn)出較小的偏差，說明高斯函數(shù)在較大的高度時擬合垂直剖面ZTD效果更佳。由表1可以看出，4個格網(wǎng)點的平均bias均為負值，說明高斯型函數(shù)擬合垂直剖面ZTD具有一定的系統(tǒng)誤差，但僅約-0.6 mm。4個格網(wǎng)點的均方根誤差分別為11.439，8.864，7.573，7.592 mm，總體來看，采用高斯函數(shù)對垂直剖面ZTD進行擬合具有較高的精度和穩(wěn)定性。

為了使式(6)能夠?qū)θ我飧叨鹊腪TD進行高程改正，將參考高程的ZTD插值到目標高程，因此，對式(6)進行推導(dǎo)[16]：

(7)

式中，e為自然常數(shù)；hr和ht分別為參考高程和目標高程，單位km；b和c為模型系數(shù)，單位km；ZTDr為參考高程hr處的ZTD值，單位m；ZTDt為目標高程ht處的ZTD值，單位m。利用式(7)即可將參考高程的ZTD插值到目標高程，更有利于實際應(yīng)用。

2 ZTD模型系數(shù)的時空特性分析

為了探究模型系數(shù)b和模型系數(shù)c的周期性，分析二者的精細時間變化特性，根據(jù)廣西地區(qū)2016—2017年ERA5大氣再分析資料積分計算的分層ZTD，采用高斯函數(shù)擬合出該研究區(qū)域內(nèi)所有格網(wǎng)點逐小時(UTC時)的模型系數(shù)b與模型系數(shù)c。選取具有代表性的一個格網(wǎng)點(25.75°N，107.5°E)，計算該格網(wǎng)點2016—2017年的日均模型系數(shù)b與模型系數(shù)c，對二者進行年周期和半年周期的三角函數(shù)擬合，并采用快速傅里葉變換分析周期性，結(jié)果如圖2所示。

由圖2(a)可知，模型系數(shù)b的值為-100～-40 km，具有顯著的季節(jié)變化特征，具體表現(xiàn)為春冬季節(jié)較大，夏季較小，且具有明顯的周期性；由圖2(c)的頻譜分析可以進一步看出，模型系數(shù)b具有顯著的年周期和半年周期。由圖2(b)可知，模型系數(shù)c的值為25～40 km，模型系數(shù)c同樣具有顯著的季節(jié)變化特征，具體表現(xiàn)為春冬季節(jié)較小，而夏季較大，同樣具有明顯的周期性；由圖2(d)的頻譜分析可以進一步看出，模型系數(shù)c也具有顯著的年周期和半年周期。

(a)b的日均時間序列及其擬合線

為了分析模型系數(shù)b與模型系數(shù)c更精細的時間變化特征，同樣選取上述格網(wǎng)點，計算年均UTC00:00—23:00的逐小時模型系數(shù)b與模型系數(shù)c，對其進行日周期和半日周期的三角函數(shù)擬合，并采用快速傅里葉變換分析周期性，結(jié)果如圖3所示。

(a)b的年均逐小時時間序列及其擬合線

由圖3(a)可知，模型系數(shù)b具有顯著的日變化特征，具體表現(xiàn)為從00:00(UTC時)開始，先增大后減小最后又增大的變化趨勢，其值在夜晚較大，白天較小,且具有明顯的周期性;由圖3(c)的頻譜分析可以進一步看出，模型系數(shù)b具有顯著的日周期和半日周期。由圖3(b)可知，模型系數(shù)c同樣具有顯著的日變化特征，具體表現(xiàn)為從00:00(UTC時)開始，先減小后增大最后又減小的變化趨勢，其值在夜晚較小，白天較大,同樣具有明顯的周期性;由圖3(d)的頻譜分析可以進一步看出，模型系數(shù)b也具有顯著的日周期和半日周期。因此，在采用高斯型函數(shù)構(gòu)建垂直剖面ZTD模型時，需要考慮模型系數(shù)b與模型系數(shù)c的年周期、半年周期、日周期以及半日周期。

為了分析模型系數(shù)b與模型系數(shù)c的年均值、年周期振幅、半年周期振幅、日周期振幅及半日周期振幅在廣西地區(qū)的區(qū)域分布特征，對每個格網(wǎng)點的模型系數(shù)b與模型系數(shù)c進行顧及年周期、半年周期、日周期及半日周期擬合，并計算每個格網(wǎng)點模型系數(shù)b與模型系數(shù)c各自的年均值、年周期振幅、半年周期振幅、日周期振幅及半日周期振幅。對離散格網(wǎng)點上的數(shù)值進行水平插值，結(jié)果如圖4和圖5所示。

由圖4可以看出顯著的區(qū)域性，其中，模型系數(shù)b年均值在廣西西北地區(qū)較大，在廣西東南地區(qū)較小，且呈現(xiàn)由西北向東南逐漸減小的變化趨勢，該變化可能與廣西西北地區(qū)海拔高、東南地區(qū)海拔低有關(guān)；年周期振幅由北向南逐漸減小，該變化可能由廣西地區(qū)的南部臨海氣候和北部內(nèi)陸氣候差異所致。而半年周期振幅、日周期振幅和半日周期振幅沒有明顯的區(qū)域變化規(guī)律，半年周期振幅在廣西西南地區(qū)較大，在東南地區(qū)較??；日周期振幅在西北和西南地區(qū)較大，在中部較小；半日周期振幅在西南地區(qū)較大，在東北地區(qū)較小。由圖5可看出，模型系數(shù)c的年均值及各振幅的區(qū)域分布特征與模型系數(shù)b的區(qū)域分布特征有極高的相似性。模型系數(shù)c的年均值與年周期振幅同樣具有顯著的區(qū)域變化規(guī)律，而半年周期振幅、日周期振幅和半日周期振幅也沒有明顯的區(qū)域變化規(guī)律。由于模型系數(shù)b與模型系數(shù)c的年均值與各振幅在不同地理位置存在較大的差異，因此，在采用高斯函數(shù)構(gòu)建垂直剖面ZTD模型時，每個格網(wǎng)點均需要考慮模型系數(shù)b與模型系數(shù)c的年周期、半年周期、日周期以及半日周期。

(a)年均值

(a)年均值

為了分析模型系數(shù)b和模型系數(shù)c與經(jīng)度的相關(guān)性，選取23°N與25°N兩個緯度上2017年的模型系數(shù)b和模型系數(shù)c，計算每個格網(wǎng)點的年均b值與年均c值，結(jié)果如圖6所示。

(a)23°N

為了分析模型系數(shù)b和模型系數(shù)c與緯度的相關(guān)性，選取108°E與111°E兩個經(jīng)度上2017年的模型系數(shù)b和模型系數(shù)c，計算每個格網(wǎng)點的年均b值與年均c值，結(jié)果如圖7所示。

(a)108°E

圖6(a)和(b)分別為模型系數(shù)b在23°N和25°N兩個緯度上隨經(jīng)度變化的散點分布及其擬合線，從圖中可看出，模型系數(shù)b與經(jīng)度呈現(xiàn)出較為顯著的負相關(guān)性；圖6(c)和(d)分別為模型系數(shù)c在23°N和25°N隨經(jīng)度變化的散點分布及其擬合線，從圖中可看出，模型系數(shù)c與經(jīng)度呈現(xiàn)出較為顯著的正相關(guān)性。圖7(a)和(b)分別為模型系數(shù)b在108°E和111°E兩個經(jīng)度上隨緯度變化的散點分布及其擬合線，從圖中可看出，模型系數(shù)b與緯度呈現(xiàn)出較為顯著的正相關(guān)性；圖7(c)和(d)分別為模型系數(shù)c在108°E和111°E兩個經(jīng)度上隨緯度變化的散點分布及其擬合線，從圖中可看出，模型系數(shù)c與緯度呈現(xiàn)出較為顯著的負相關(guān)性。對比圖6和圖7可知，模型系數(shù)b、模型系數(shù)c與緯度的相關(guān)性更高。

為了進一步分析模型系數(shù)b、模型系數(shù)c與經(jīng)緯度的相關(guān)性，利用皮爾遜相關(guān)性分析法計算b值、c值與經(jīng)緯度的相關(guān)系數(shù)，結(jié)果如表2所示。

表2 模型系數(shù)b/c與經(jīng)緯度的相關(guān)系數(shù)

在皮爾遜相關(guān)性分析法中，2組數(shù)據(jù)的相關(guān)系數(shù)越接近1或-1，說明2組數(shù)據(jù)相關(guān)性越高，相關(guān)系數(shù)為正時，說明2組數(shù)據(jù)呈正相關(guān)，相關(guān)系數(shù)為負時，說明2組數(shù)據(jù)呈負相關(guān)。由表2可知，模型系數(shù)b與經(jīng)度的相關(guān)系數(shù)為-0.824，二者呈負相關(guān)關(guān)系，相關(guān)性較高；模型系數(shù)b與緯度的相關(guān)系數(shù)為0.969，二者呈正相關(guān)關(guān)系，相關(guān)性較高，相比經(jīng)度，模型系數(shù)b與緯度的相關(guān)性更高。模型系數(shù)c與經(jīng)度的相關(guān)系數(shù)為0.822，二者呈正相關(guān)關(guān)系，相關(guān)性較高；模型系數(shù)c與緯度的相關(guān)系數(shù)為-0.968，二者呈負相關(guān)關(guān)系，相關(guān)性較高，相比經(jīng)度，同樣可以看出模型系數(shù)c與緯度的相關(guān)性更高。由于模型系數(shù)b、模型系數(shù)c與經(jīng)緯度均存在較高的相關(guān)性，因此，在采用高斯型函數(shù)構(gòu)建垂直剖面ZTD模型時，均需要考慮模型系數(shù)b、模型系數(shù)c與經(jīng)緯度的相關(guān)性。

3 結(jié)束語

針對廣西地區(qū)地形起伏較大、氣候復(fù)雜的特點，已有的低時空分辨率ZTD垂直剖面模型難以滿足該地區(qū)的需求，根據(jù)2016—2017年廣西地區(qū)ERA5大氣再分析資料，分析基于高斯函數(shù)的ZTD垂直剖面模型的模型系數(shù)b與模型系數(shù)c的時空分布特征，結(jié)果表明：

① 模型系數(shù)b與模型系數(shù)c均存在季節(jié)變化特征，二者季節(jié)變化差異較大，模型系數(shù)b的數(shù)值表現(xiàn)出春冬季節(jié)較大而夏季較小的特征，模型系數(shù)c的數(shù)值表現(xiàn)出春冬季節(jié)較小而夏季較大的特征，二者均存在顯著的年周期與半年周期。

② 模型系數(shù)b與模型系數(shù)c在日變化方面表現(xiàn)出一定的差異性，模型系數(shù)b的數(shù)值在夜晚較大而在白天較小，模型系數(shù)c的數(shù)值與模型系數(shù)b的情況相反，二者均存在顯著的日周期和半日周期。

③ 模型系數(shù)b與模型系數(shù)c的年均值、年周期振幅、半年周期振幅、日周期振幅及半日周期振幅在廣西地區(qū)存在顯著的區(qū)域性特征，其中，b值、c值的年均值與年周期振幅存在沿某一方向逐漸變化的規(guī)律。

④ 模型系數(shù)b、模型系數(shù)c與經(jīng)緯度均存在較高的相關(guān)性，二者與緯度的相關(guān)性高于與經(jīng)度的相關(guān)性。

綜上可知，在廣西地區(qū)構(gòu)建基于高斯函數(shù)的ZTD垂直剖面模型時需考慮模型系數(shù)b與模型系數(shù)c的年周期、半年周期、日周期、半日周期，以及二者與經(jīng)緯度的相關(guān)性。研究結(jié)果可為構(gòu)建基于高斯函數(shù)的ZTD垂直剖面模型提供參考。