江俊飛，王磊，葉世榕，羅穎婷，李大江，鄂盛龍，賴詩鈺

(1.廣東電網(wǎng)有限責任公司電力科學研究院，廣東 廣州 510080；2.國家衛(wèi)星定位系統(tǒng)工程技術研究中心，湖北 武漢 430000；3.中山大學 測繪科學與技術學院，廣東 珠海 519080)

目前，無人機自動巡檢技術已成為電網(wǎng)公司開展設備巡視和運維的主要手段之一[1]，其依賴于全球衛(wèi)星導航系統(tǒng)(global navigation satellite system，GNSS)和地面服務系統(tǒng)。為了能保障無人機定位的穩(wěn)定性和可靠性，廣東電網(wǎng)有限責任公司自建了連續(xù)運行參考站(continuous operating reference station，CORS)定位服務系統(tǒng)；但在實際應用中受廣東地區(qū)地理位置影響，無人機高精度定位仍存在一定的定位誤差。其中，對流層天頂總延遲(zenith total delay，ZTD)是定位大氣誤差中的主要誤差，影響定位精度可達至分米級別。ZTD產(chǎn)生的原因是衛(wèi)星信號在穿越對流層過程中受水汽等大氣介質影響，傳播路徑和傳播速度發(fā)生改變。通常情況下，ZTD分為靜水延遲和濕延遲，濕延遲作為受大氣水汽影響較大的部分，是ZTD精準估算的難點。采用常規(guī)的Hopfiled模型[2]與Saastamonien模型[3]獲取ZTD，在水汽影響較大區(qū)域其精度較差，并且依賴于氣象站實測氣象數(shù)據(jù)，成本較高。為避免采用實測氣象數(shù)據(jù)，有學者通過插值表或球諧函數(shù)構建氣象參數(shù)模型，如Boehm等人的GPT系列的模型[4-6]、Collins等人的UNB3[7]和UNB3m[8]模型、Penna等人的EGONS模型[9]，以及Mateus等人通過對ERA5再分析資料進行時序分析而最新提出的HGPT模型[10]和HGPT2模型[11]。為了更加直接地構建ZTD模型，并且避免通過參數(shù)計算ZTD引入的偏差，部分學者直接針對經(jīng)高精度解算或數(shù)值氣象資料解算的長期ZTD周期性變化來構建模型，如Krueger等構建的TropGrid模型[12]，Schüler等人改進的球諧函數(shù)模型TropGrid2[13]，Yao等人構建的ITG模型[14]、GZTD模型[15]、GZTD2模型[16]，以及黃良珂等人構建的GGZTD模型[17]。上述模型主要基于歷史ZTD的時間序列，能反映出對流層延遲的長期變化趨勢，誤差為4～6 cm；但其很難反映對流層延遲的短期變化，且在不同區(qū)域的表現(xiàn)效果并不相同，在一些特定的、濕延遲受水汽影響較大的地區(qū)，這些經(jīng)驗模型的精度影響實時精密單點定位的快速解算以及ZTD反演水汽的準確度。同時，近些年厄爾尼諾現(xiàn)象的頻繁發(fā)生，較大影響了我國的沿海地區(qū)，導致水汽劇烈變化[18]。構建區(qū)域實時ZTD模型不僅對實時定位有著重要意義，也能為災害氣候預警提供一項可靠的數(shù)據(jù)源[19]。

為解決廣東地區(qū)常規(guī)ZTD計算模型精度不夠高的問題，本研究首先分析ZTD隨高程變化常用的擬合函數(shù)(指數(shù)函數(shù)、高斯函數(shù))，提出新的分段擬合方式(在0～2 km采用線性擬合，2～8 km和8～50 km分別采用指數(shù)函數(shù)擬合)，并以廣東地區(qū)2020—2022年的ERA5數(shù)據(jù)和2023年1—3月的GNSS站點數(shù)據(jù)實時解算ZTD，構建修正后的區(qū)域實時ZTD模型；接著，根據(jù)3年ERA5數(shù)據(jù)，使用該模型計算得到分段擬合函數(shù)系數(shù)的時間序列，并利用計及年、半年和日周期變化的傅里葉級數(shù)擬合該時間序列，獲得ZTD經(jīng)驗格網(wǎng)模型；然后，計算用經(jīng)驗模型在GNSS站點位置處求出的ZTD與站點實時GNSS-ZTD的偏差，利用反距離加權的方法計算出經(jīng)驗格網(wǎng)點的實時偏差，并補償格網(wǎng)點處的經(jīng)驗擬合函數(shù)，獲得具有實時特性的ZTD格網(wǎng)數(shù)據(jù)；最后，利用廣東電網(wǎng)CORS服務系統(tǒng)歷史數(shù)據(jù)，對獲得的格網(wǎng)數(shù)據(jù)進行精度驗證。

1 基于ERA5再分析資料的廣東地區(qū)ZTD精細建模

ERA5是由歐洲中期天氣預報中心(European Centre for Medium-range Weather Forecasts，ECMWF)所研發(fā)和維護的全球氣象重分析數(shù)據(jù)集，以高分辨率和豐富的氣象參數(shù)描述地球大氣的時空變化[10]。該數(shù)據(jù)集的空間分辨率為0.25°×0.25°，時間分辨率為1 h，其包含了自1979年以來的溫度、濕度、風速、氣壓、降水等多種氣象參數(shù)，不僅涵蓋了地球表面的情況，還考慮了不同垂直層次的大氣特性。ERA5數(shù)據(jù)以NetCDF格式提供。本文采用的ERA5數(shù)據(jù)空間分辨率為0.5°×0.5°，共計37層氣壓分層數(shù)據(jù)，時間分辨率為1 h(下載地址：https：//cds.climate. copernicus.eu/)。首先根據(jù)ERA5產(chǎn)品中氣象參數(shù)(氣壓、水汽壓、溫度)采用方程計算出每層的大氣折射率，再從各層級所處的位置對折射率進行積分求和，最后加上ERA5最高層以上的ZTD，便可以得到每個氣壓層對應的ZTD(量符號記為KZTD)：

(1)

(2)

式(1)、(2)中：N為大氣折射率；e為為水汽壓，單位hpa；T為溫度，單位K；p為氣壓，單位hpa；φ為格網(wǎng)點所在緯度；H為不同氣壓層的大地高，H1、Hn分別為第1、n層氣壓層的大地高，單位km。將頂層氣壓p代入Saastamonien經(jīng)驗模型，計算得到ERA5氣象數(shù)據(jù)頂層以上的ZTD，因50 km以上大氣的水汽含量基本為0，經(jīng)驗公式能夠精確計算該高度上的ZTD。經(jīng)過計算后得到的ERA5數(shù)據(jù)字段包括時間、經(jīng)緯度、37層氣壓層的大地高，以及每層的ZTD。

實時站點數(shù)據(jù)是使用gamp實時計算出的ZTD數(shù)據(jù)，數(shù)據(jù)包括時間、經(jīng)緯度、大地高和ZTD。

GPT3[6]在ECMWF建立的模型基礎上增加了梯度模型，改良了映射函數(shù)關系，是目前主流的先驗對流層延遲模型之一。該模型只需提供用戶坐標和年積日即可獲得所需要的氣象參數(shù)，最后利用Saastamoinen模型獲得用戶處的ZTD。

HGPT2[11]同樣也是一種經(jīng)驗模型，主要是對全球溫度、氣壓和相對濕度的時間序列進行建模，通過自身擬合的氣象參數(shù)和Saastamoinen模型，來獲得全球任意時空的ZTD數(shù)據(jù)。

1.1 ZTD垂直剖面模型高程歸化擬合函數(shù)的改進

目前，ZTD常用的垂直擬合函數(shù)有指數(shù)函數(shù)

KZTD=aexpbH

(3)

和高斯函數(shù)[20-21]

KZTD=aexp[-(bH)2/c2].

(4)

式(3)、(4)中a、b、c均為擬合系數(shù)。隨機選取不同經(jīng)緯度上的格網(wǎng)點，分別采用指數(shù)函數(shù)和高斯函數(shù)進行非線性最小二乘擬合，擬合結果如圖1所示。

圖1 在不同高度的指數(shù)函數(shù)和高斯函數(shù)擬合ZTDFig.1 Exponential and Gaussian functions fit ZTD at different heights

為了進一步觀察這2種擬合函數(shù)在不同高度處的擬合效果，統(tǒng)計經(jīng)度73°～133°、緯度3°～53°所有格網(wǎng)點不同高度處的ZTD擬合偏差，其0～50 km區(qū)間的偏差曲線如圖2所示，將ZTD的擬合偏差記為“DZTD”。

圖2 在不同高度的指數(shù)函數(shù)和高斯函數(shù)擬合偏差Fig.2 Fitting bias of exponential and Gaussian functions at different heights

從圖2可以看出：指數(shù)函數(shù)在0～<2 km、2～<8 km、10～<30 km均會出現(xiàn)較大的偏差，最大偏差可達到6 cm；高斯函數(shù)在0～6 km會出現(xiàn)較大的偏差，最大也可達到6 cm，但其在高空區(qū)域的擬合效果較好，基本都在2 cm以內。無論采用指數(shù)函數(shù)還是高斯函數(shù)，在擬合ZTD與高程的變化關系時，都會在較低高層范圍內出現(xiàn)較大的偏差。

通常情況下，GNSS觀測站測量高程在3 km以下。經(jīng)過大量實驗，本研究采用如下分段擬合函數(shù)(以下簡稱“分段函數(shù)”)來表示不同高度區(qū)間內ZTD的變化：

(5)

式(5)中a1—a3以及b1—b3均為擬合形式中的未知參數(shù)。從圖1可以看出，H在0～<2 km時，ZTD隨高程基本是線性變化，故本研究在0～<2 km時使用線性函數(shù)；同時由圖2可以看出，指數(shù)函數(shù)擬合的偏差在2～<8 km和20 km左右較大，故在2～<8 km和8 km以上使用指數(shù)函數(shù)分別進行非線性最小二乘擬合。在經(jīng)度85°、緯度40°位置，分段函數(shù)擬合效果和偏差統(tǒng)計如圖3所示。

圖3 在不同高度的分段函數(shù)的擬合ZTD及DZTDFig.3 Segmented function fits ZTD and DZTD at different heights

從圖3可以看出，新的分段函數(shù)效果較采用指數(shù)函數(shù)/高斯函數(shù)更好，偏差基本在1 cm左右。分段函數(shù)的擬合形式避免了高斯函數(shù)和指數(shù)函數(shù)在低空出現(xiàn)較大偏差的情況。統(tǒng)計3種擬合函數(shù)在不同高度范圍擬合ZTD的均方根誤差(root mean square error，RMSE)，通過RMSE來表征精度，見表1。

表1 在不同高度范圍的指數(shù)函數(shù)、高斯函數(shù)、分段函數(shù)的RMSETab.1 RMSE values of exponential，Gaussian，and segmental fitting functions in different height ranges

在0～<5 km的高度范圍內，分段函數(shù)展現(xiàn)出最低的RMSE，為0.6 cm，相比之下，指數(shù)函數(shù)、高斯函數(shù)的RMSE分別為1.1 cm、1.5 cm；類似地，在5～<10 km、10～<20 km以及20～<50 km的高度范圍內，分段函數(shù)也呈現(xiàn)出最低的RMSE，分別為0.5 cm、0.7 cm和0.3 cm。總體而言，在任意高度范圍內，分段函數(shù)都表現(xiàn)出較低的誤差。分段函數(shù)的0～50 km全高度RMSE(簡稱為“總RMSE”)為0.6 cm，而高斯函數(shù)、指數(shù)函數(shù)的總RMSE分別為1.1 cm、1.4 cm。

以東經(jīng)105°和北緯35°為分界線將中國區(qū)域劃分成東北地區(qū)、東南地區(qū)、西北地區(qū)、西南地區(qū)，分別統(tǒng)計這3種擬合函數(shù)在不同地區(qū)的總RMSE，見表2。

表2 3種擬合函數(shù)在不同地區(qū)的總RMSETab.2 RMSE values of three fitted functions in different regions

由表2可看出：在東北地區(qū)，分段函數(shù)表現(xiàn)較好，其總RMSE達到0.5 cm，指數(shù)函數(shù)、高斯函數(shù)的總RMSE分別為1.2 cm、1.1 cm；在東南地區(qū)，分段函數(shù)的總RMSE為0.6 cm，高斯函數(shù)的總RMSE(1.4 cm)略大于指數(shù)函數(shù)的總RMSE(1.0 cm)；西北地區(qū)中，高斯函數(shù)和分段函數(shù)呈現(xiàn)出相同的精度，總RMSE均為0.4 cm，指數(shù)函數(shù)的總RMSE為1.9 cm，指數(shù)函數(shù)在西北地區(qū)的擬合精度較低；至于西南地區(qū)，指數(shù)函數(shù)、高斯函數(shù)顯示出接近的精度，總RMSE分別為1.5 cm、1.2 cm，分段函數(shù)的總RMSE為0.6 cm。由此可知：分段函數(shù)在不同區(qū)域的擬合效果都比較穩(wěn)定，且總RMSE一直保持在0.5 cm左右；指數(shù)函數(shù)在東北、東南地區(qū)的精度比西北、西南地區(qū)的精度要較高一些；高斯函數(shù)在西北地區(qū)精度較高，其余3個地區(qū)的精度要稍差一些。

1.2 構建廣東地區(qū)ZTD格網(wǎng)經(jīng)驗模型

廣東地區(qū)位于東南沿海地區(qū)，根據(jù)上述分析結果，采用分段函數(shù)來對該地區(qū)進行建模，擬合函數(shù)見式(5)。使用積分法計算出廣東地區(qū)ERA5數(shù)據(jù)2020年1月—2022年1月的、每個格網(wǎng)點不同氣壓層的ZTD數(shù)據(jù)，再使用式(5)分段函數(shù)計算出每個時刻的擬合函數(shù)系數(shù)，得到3年擬合函數(shù)系數(shù)的時間序列。圖4所示為格網(wǎng)點(緯度23°、經(jīng)度113°)擬合系數(shù)a1—a3、b1—b3的時間序列，曲線中各數(shù)據(jù)點對應各年份中整點時刻的數(shù)據(jù)。

圖4 2020年至2022年分段函數(shù)系數(shù)的時間序列Fig.4 Time series of coefficients of segmented fitting functions from 2020 to 2022

從圖4可以看出這些系數(shù)具有較強的周期性變化。以系數(shù)a1為例，對系數(shù)的時間序列進行傅里葉頻譜分析(其中采樣頻率設為1 Hz)，結果如圖5所示。

圖5 系數(shù)a1的頻譜分析Fig.5 Spectral analysis of coefficient a1

由圖5得到擬合系數(shù)的時間序列主要頻率，通過傅里葉頻譜分析，最終對每個系數(shù)采用計及年、半年和日周期變化的傅里葉級數(shù)，擬合計算式為

(6)

式中：k、ki(i=1，…，6)均為模型系數(shù)，k表示均值項，k1、k2表示年周期項，k3、k4表示半年周期項，k5、k6表示日周期項；tdoy為年積日；thod為日內時；m為分段函數(shù)的系數(shù)，m∈{a1，a2，a3，b1，b2，b3}。使用式(6)和非線性最小二乘方法，擬合并保存廣東地區(qū)每個格網(wǎng)點擬合系數(shù)的3年時間序列數(shù)據(jù)，構建完成廣東地區(qū)ZTD經(jīng)驗模型。輸入時間(tdoy和thod)便可以得到任意格網(wǎng)點某個時刻不同高度處的ZTD數(shù)據(jù)。該模型的空間分辨率為0.5°×0.5°，時間分辨率1 h，區(qū)域范圍為緯度20°～26°、經(jīng)度109°～118°。

1.3 實時解算GNSS修正ERA5經(jīng)驗模型

經(jīng)驗模型能反映ZTD的長期變化趨勢，但很難反映其短期實時變化，因此可以加入實時站點數(shù)據(jù)。步驟如下：

a)輸入當前時刻的年積日和日內時，獲得當前時刻每個格網(wǎng)點的ZTD隨高程變化的擬合函數(shù)；

b)利用雙線性插值得到GNSS站點處的經(jīng)驗ZTD，并與GNSS站的實時ZTD作差，得到經(jīng)驗模型的實時DZTD；

c)搜尋每個格網(wǎng)點距離最近的3個GNSS站點，計算出站點DZTD，然后采用反距離加權的方法獲得格網(wǎng)點處的DZTD(將該DZTD記為“DZTDg”)；

d)將格網(wǎng)點處的擬合函數(shù)加上實時DZTDg，便能得到實時的ZTD格網(wǎng)(將該格網(wǎng)模型記為“Gdgrid”)數(shù)據(jù)，其格式為時間、經(jīng)緯度、ZTD隨高程的擬合系數(shù)，以及計算出來的DZTD。

用戶使用步驟：

a)接收到播發(fā)的Gdgrid數(shù)據(jù)后，找到所處格網(wǎng)；

b)輸入高程獲得4個格網(wǎng)點在用戶高程處的ZTD；

c)利用雙線性插值計算得到用戶所在位置ZTD。

2 模型精度驗證

按照實時精密單點定位數(shù)據(jù)處理方法，計算獲得2023年1—3月廣東電網(wǎng)96個GNSS站的實時ZTD數(shù)據(jù)。計算處理策略見表3。

表3 精密單點定位參數(shù)設置Tab.3 Parameter setting of PPP

選用廣東電網(wǎng)31個GNSS站作為基準站，65個GNSS站作為驗證站。選取的基準站和驗證站分布均勻，能較好反映不同地區(qū)的精度差異。對于參與建模的站，按上述方法對每個格網(wǎng)點進行操作得到Gdgrid數(shù)據(jù)，并利用驗證站對Gdgrid數(shù)據(jù)精度進行評定，部分GNSS驗證站的總RMSE見表4。

從表4可以看出，在這25個驗證站上Gdgrid的總RMSE都在1 cm以下，僅有個別站總RMSE在1 cm以上。通過觀察驗證站位置后發(fā)現(xiàn)，導致這些站精度較差的原因是這些驗證站多數(shù)位于邊緣地帶，僅有一側具有基準站，或者是距離基準站較遠。GPT3、HGPT2在部分站精度較高，但在有些站會有較大的偏差，原因是GPT3和HGPT2屬于經(jīng)驗模型，并沒有添加實時數(shù)據(jù)源，導致在部分地區(qū)這2個模型的精度較差。統(tǒng)計在所有驗證站Gdgrid的總RMSE分布情況，如圖6所示。

圖6 驗證站2023年1—3月Gdgrid的總RMSEFig.6 RMSE values of Gdgrid of validation stations from January to March 2023

由圖6可以看出，僅有5個站Gdgrid的總RMSE在1.5 cm以上，多數(shù)站總RMSE都在1.0 cm以下，可見本文方法獲得的實時ZTD數(shù)據(jù)精度較高。

統(tǒng)計所有驗證站在3個月內的總RMSE，并與GPT3、HGPT2經(jīng)驗模型的總RMSE進行比較，見表5。

表5 Gdgrid與GPT3、HGPT2的對比Tab.5 Comparison of GPT3，HGPT2 and Gdgrid

從表5可以看出Gdgrid的總RMSE僅有1.12 cm，遠優(yōu)于GPT3、HGPT2的總RMSE，其中HGPT2精度最差。統(tǒng)計結果表明，相對于經(jīng)驗模型，附加實時數(shù)據(jù)的修正ERA5的模型精度更高，能夠更好反映ZTD的短期變化。

3 結論

本研究針對現(xiàn)有基于指數(shù)函數(shù)、高斯函數(shù)的ZTD模型在低空范圍會出現(xiàn)較大偏差的問題，提出基于分段擬合的新建模方法，并采用實際觀測數(shù)據(jù)對模型精度進行了驗證，結論如下：

a)本文提出的ZTD模型在海拔高度0～<2 km采用線性函數(shù)，2～<8 km和8 km以上分別采用不同擬合系數(shù)的指數(shù)函數(shù)。結果表明分段函數(shù)的擬合精度最高，尤其是在0～10 km高度范圍，指數(shù)函數(shù)和高斯函數(shù)的RMSE基本在1 cm以上，而分段函數(shù)的RMSE在0.5 cm左右，極大提高了ZTD在低空區(qū)域的擬合精度。

b)實驗結果表明：分段函數(shù)在中國不同區(qū)域處的擬合精度基本不變，總RMSE維持在0.5 cm左右；指數(shù)函數(shù)在西北地區(qū)精度最差，總RMSE為1.9 cm，在東北、東南、西南地區(qū)的總RMSE在1.0～1.5 cm之間；高斯函數(shù)在西北地區(qū)表現(xiàn)較好，總RMSE僅有0.4 cm，在東北、東南、西南地區(qū)精度與指數(shù)函數(shù)相差不大。

c)根據(jù)廣東電網(wǎng)96個GNSS站2023年1—3月的實時ZTD數(shù)據(jù)，對Gdgrid模型格網(wǎng)數(shù)據(jù)進行精度驗證，結果表明僅有5個站總RMSE大于1.5 cm，其余站的總RMSE基本在1.0 cm以下。

d)對比Gdgrid與GPT3、HGPT2模型精度，Gdgrid的精度遠優(yōu)于GPT3、HGPT2，Gdgrid的總RMSE為1.12 cm，而GPT3、HGPT2的總RMSE分別為4.24 cm和5.13 cm。