基于高斯函數的ZTD模型系數時空特性分析

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

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

0 引言

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

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

1 數據來源與計算原理

1.1 數據來源

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

1.2 計算原理

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

(1)

(2)

(3)

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

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

(4)

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

(5)

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

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

(6)

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

(a)26°N,111°E

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

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

為了使式(6)能夠對任意高度的ZTD進行高程改正，將參考高程的ZTD插值到目標高程，因此，對式(6)進行推導[16]：

(7)

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

2 ZTD模型系數的時空特性分析

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

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

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

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

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

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

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

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

(a)年均值

(a)年均值

為了分析模型系數b和模型系數c與經度的相關性，選取23°N與25°N兩個緯度上2017年的模型系數b和模型系數c，計算每個格網點的年均b值與年均c值，結果如圖6所示。

(a)23°N

為了分析模型系數b和模型系數c與緯度的相關性，選取108°E與111°E兩個經度上2017年的模型系數b和模型系數c，計算每個格網點的年均b值與年均c值，結果如圖7所示。

(a)108°E

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

為了進一步分析模型系數b、模型系數c與經緯度的相關性，利用皮爾遜相關性分析法計算b值、c值與經緯度的相關系數，結果如表2所示。

表2 模型系數b/c與經緯度的相關系數

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

3 結束語

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

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

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

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

④ 模型系數b、模型系數c與經緯度均存在較高的相關性，二者與緯度的相關性高于與經度的相關性。

綜上可知，在廣西地區構建基于高斯函數的ZTD垂直剖面模型時需考慮模型系數b與模型系數c的年周期、半年周期、日周期、半日周期，以及二者與經緯度的相關性。研究結果可為構建基于高斯函數的ZTD垂直剖面模型提供參考。