張 爽，陳西宏，劉 強(qiáng)，劉 贊,3，王慶力

1. 空軍工程大學(xué)防空反導(dǎo)學(xué)院，陜西 西安 710051； 2. 93305部隊(duì)，遼寧 沈陽(yáng) 110000； 3. 93567部隊(duì)，河北 保定 074100

精準(zhǔn)的對(duì)流層天頂延遲估計(jì)對(duì)導(dǎo)航定位授時(shí)和甚長(zhǎng)基線干涉測(cè)量有重要影響[1-3]?；贜WM和折射率積分法可獲取全球任意位置ZTD，具有計(jì)算精度高、可用范圍廣、不受時(shí)空限制等優(yōu)點(diǎn)，是目前分析ZTD時(shí)空分布、構(gòu)建和驗(yàn)證ZTD經(jīng)驗(yàn)?zāi)Ｐ偷闹匾侄沃籟4-10]。同時(shí)，隨著全球氣象同化觀測(cè)資料的不斷豐富和同化手段的進(jìn)步，NWM時(shí)空分辨率和數(shù)據(jù)質(zhì)量在不斷提高。以2019年歐洲中期天氣預(yù)報(bào)中心(European Centre For Medium Range Weather Forecasts，ECMWF)發(fā)布的ERA5氣壓分層產(chǎn)品為例，其空間分辨率可達(dá)0.25°×0.25°，時(shí)間分辨率為1 h，可提供高質(zhì)量的氣象數(shù)據(jù)。NWM將全球氣象觀測(cè)數(shù)據(jù)通過物理模型基于數(shù)據(jù)同化原理生成[11]，其數(shù)據(jù)精度優(yōu)于氣象經(jīng)驗(yàn)?zāi)Ｐ停源嬖谝欢ㄕ`差。目前常用的精度評(píng)估方法是利用GNSS ZTD產(chǎn)品進(jìn)行外部評(píng)估。文獻(xiàn)[8，12—15]利用全球或區(qū)域性GNSS測(cè)站ZTD產(chǎn)品對(duì)諸如ECMWF提供的ERA-interim和ERA5氣壓分層產(chǎn)品、美國(guó)國(guó)家大氣研究中心(National Center for Atmospheric Research，NCAR)的WRF中尺度數(shù)值氣象預(yù)報(bào)模型、美國(guó)國(guó)家海洋和大氣管理局(National Oceanic and Atmospheric Administra-tion，NOAA)的GDAS數(shù)值氣象模型進(jìn)行精度評(píng)估。這種評(píng)估方法依賴于GNSS站點(diǎn)分布，只能給出測(cè)站點(diǎn)處的ZTD精度，對(duì)GNSS站點(diǎn)未配置區(qū)域只能以區(qū)域性平均值代替，誤差較大。

考慮到NWM誤差主要來源于同化觀測(cè)資料的空間分布、同化觀測(cè)資料數(shù)據(jù)的時(shí)間演變及局部天氣變化引起的不確定性[16]，而局部天氣變化一方面可通過NWM自身氣象數(shù)據(jù)部分反應(yīng)，另一方面又與地理環(huán)境特征存在相關(guān)性，因此本文以ECMWF發(fā)布的ERA5氣壓分層產(chǎn)品為研究對(duì)象，通過NWM自身氣象數(shù)據(jù)和地形特征數(shù)據(jù)構(gòu)建特征向量，并基于粒子群算法和擴(kuò)展的RBF神經(jīng)網(wǎng)絡(luò)構(gòu)建ZTD精度估計(jì)模型，再以我國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)(Crustal Movement Observation Network of China,CMONOC，簡(jiǎn)稱“陸態(tài)網(wǎng)絡(luò)”)的GNSS測(cè)站ZTD產(chǎn)品作為參考獲取ZTD參考精度，結(jié)合特征向量訓(xùn)練網(wǎng)絡(luò)模型，實(shí)現(xiàn)不依賴外部數(shù)據(jù)的ZTD精度估計(jì)。

1 數(shù)據(jù)集生成與劃分

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

數(shù)值氣象模型采用ECMWF提供的ERA5氣壓分層產(chǎn)品(空間分辨率為0.25°、時(shí)間分辨率為1 h)，截取2016—2020年我國(guó)及其周邊區(qū)域(70°E—140°E，0°N—55°N)氣壓、溫度、比濕度和位勢(shì)作為基礎(chǔ)數(shù)據(jù)。ZTD參考值使用中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)(http:∥www.cgps.ac.cn)提供的CMONOC的GNSS測(cè)站對(duì)流層天頂延遲產(chǎn)品(時(shí)間分辨率為1 h[17]，部分測(cè)站精度經(jīng)驗(yàn)證可達(dá)4 mm[14])，從260個(gè)GNSS基準(zhǔn)站中選取月數(shù)據(jù)完整度高于60%的246個(gè)站點(diǎn)開展模型訓(xùn)練和算法驗(yàn)證，測(cè)站分布如圖1所示。本文用于構(gòu)建地形特征的高程數(shù)據(jù)來源于NASA提供的SRTM3全球數(shù)字高程數(shù)據(jù)模型經(jīng)重采樣生成的分辨率為0.25°×0.25°的數(shù)字高程地圖。

圖1 246個(gè)GNSS測(cè)站分布與海拔高度Fig.1 Distribution and altitude of 246 GNSS stations

1.2 氣象特征集

ZTD可分天頂靜力學(xué)延遲(zenith hydrostatic delay，ZHD)和天頂濕延遲(zenith wet delay，ZWD)兩部分分別計(jì)算。ZHD是氣壓和溫度的函數(shù)，現(xiàn)有研究表明無論是數(shù)值氣象模型還是基于其構(gòu)建的經(jīng)驗(yàn)?zāi)Ｐ停琙HD精度均可達(dá)到較高水平[10,18]。ZWD是溫度和水汽壓的函數(shù)，其中水汽壓分布隨機(jī)性強(qiáng)，難以準(zhǔn)確捕獲，是目前公認(rèn)的影響ZTD計(jì)算精度的主要因素。因此氣象特征集選擇地表溫度、地表水汽、ZTD月均值和ZTD月標(biāo)準(zhǔn)差作為特征向量，同時(shí)考慮到大部分地區(qū)ZTD精度存在季節(jié)性周期變化，因此將月份也作為氣象特征向量之一。

本文使用的ERA5氣壓層產(chǎn)品提供37層等壓面氣象數(shù)據(jù)，不直接提供地表數(shù)據(jù)，對(duì)地表氣壓、溫度和水汽壓采用插值或外推方式獲取。當(dāng)?shù)乇砦恢迷诘葔好孀畹讓右陨蠒r(shí)，水平方向均采用雙線性插值，垂直方向上溫度和比濕度采用線性插值，氣壓使用指數(shù)模型插值，公式為[6]

(1)

式中，P為目標(biāo)點(diǎn)處氣壓，單位為hPa；P0、T0、q0分別為鄰近網(wǎng)格點(diǎn)處氣壓、溫度(K)和比濕度；Tv為虛溫；g0為重力加速度標(biāo)準(zhǔn)值(9.806 65 m/s2)；dh為目標(biāo)點(diǎn)與網(wǎng)格點(diǎn)高程差；dMtr為干空氣摩爾質(zhì)量(28.965×10-3kg/mol)；Rd為氣體常數(shù)(8.314 3 J/K·mol)。當(dāng)?shù)乇砦恢迷诘葔好孀畹讓右韵聲r(shí)，此時(shí)垂直方向無法直接插值計(jì)算，為獲取相關(guān)參數(shù)信息，溫度值計(jì)算借鑒構(gòu)建全球氣象模型的通用方法，假設(shè)溫度以恒定溫度梯度-0.006 5℃/m[19]線性變化；比濕度則根據(jù)ERA5使用手冊(cè)推薦方法認(rèn)定其在最底層等壓面以下保持恒定不變[11]；氣壓仍利用式(1)通過指數(shù)模型計(jì)算。上述計(jì)算方法獲取的底層等壓面以下的氣象數(shù)據(jù)使用了平均經(jīng)驗(yàn)值和經(jīng)驗(yàn)?zāi)Ｐ?，算法?jiǎn)單易用，但計(jì)算精度較插值法低，考慮到一般在地表海拔低于海平面高度時(shí)會(huì)出現(xiàn)需要計(jì)算最底層等壓面以下氣象數(shù)據(jù)的情況，因此在低海拔地區(qū)需注意該方法可能帶來的額外計(jì)算誤差。

ZTD采用折射率積分法與Saastamoinen模型相結(jié)合方式，沿天頂方向分層積分至等壓面最頂層，頂層以上至對(duì)流層頂(海拔86 km)采用Saastamoinen模型計(jì)算

(2)

(3)

式中，P、q和T分別為對(duì)應(yīng)位置處的氣壓(單位為hPa)、比濕度和溫度(單位為K)；k1、k2、k3為大氣折射系數(shù)(77.60 K/hPa、69.40 K/hPa、370 100 K2/hPa)。

1.3 地形特征集

地貌地形對(duì)氣候影響體現(xiàn)在兩方面：一是占據(jù)巨大空間范圍或特殊的空間分布的大尺度宏觀地貌，如山脈、高原、盆地等，使區(qū)域內(nèi)形成獨(dú)特的氣候類型；二是中小尺度地形，如海拔高度、地面起伏度等，氣象要素分布隨地形而變化，形成局部復(fù)雜氣候。NWM采用網(wǎng)格化方式存儲(chǔ)氣象要素，在復(fù)雜或特殊地形條件下，其數(shù)據(jù)質(zhì)量勢(shì)必受到挑戰(zhàn)，因此須從局部和宏觀兩個(gè)層面構(gòu)建地形特征向量。

為確定局部地形和宏觀地貌最佳量化范圍，以高程標(biāo)準(zhǔn)差為例，分析不同范圍下高程標(biāo)準(zhǔn)差與ZTD精度相關(guān)系數(shù)。高程標(biāo)準(zhǔn)差描述了地形起伏變化情況，采用矩形窗口法計(jì)算一定區(qū)域內(nèi)的高程標(biāo)準(zhǔn)差，即獲取以目標(biāo)點(diǎn)為中心的等邊矩形所覆蓋區(qū)域地形高程標(biāo)準(zhǔn)差，并計(jì)算其與ZTD精度相關(guān)性。圖2為矩形窗口邊長(zhǎng)1000 km以下和1000 km以上時(shí)高程標(biāo)準(zhǔn)差與ZTD年均RMS相關(guān)系數(shù)。ZTD年均RMS由式(4)獲取

圖2 不同矩形窗口邊長(zhǎng)的高程標(biāo)準(zhǔn)差與ZTD年均RMS相關(guān)系數(shù)Fig.2 Correlation coefficient between elevation standard deviation of different rectangular window side lengths and ZTD mean annual RMS

(4)

由圖2可知，隨窗口面積增加，相關(guān)性呈現(xiàn)出先升高后降低，而后再升高再降低的趨勢(shì)，印證了局部地形和宏觀地貌對(duì)ZTD計(jì)算精度的多重影響，且局部地形在矩形邊長(zhǎng)230 km、宏觀地形在矩形邊長(zhǎng)4700 km時(shí)相關(guān)性最高，因此本文以矩形邊長(zhǎng)230 km和470 km作為構(gòu)建地形特征向量的基礎(chǔ)。

地形地貌難以通過單一指標(biāo)準(zhǔn)確度量，因此將高程標(biāo)準(zhǔn)差、地形切割深度、地形起伏度、地表粗糙度和高程變異系數(shù)的230 km和4700 km矩形窗口值，以及高程、經(jīng)度等多個(gè)指標(biāo)綜合起來構(gòu)建地形特征集。

1.4 輸出目標(biāo)集

模型的輸出目標(biāo)是基于ERA5氣壓分層產(chǎn)品計(jì)算的ZTD月精度，即以CMONOC的GNSS對(duì)流層天頂延遲產(chǎn)品為參考基準(zhǔn)評(píng)估在GNSS測(cè)站位置處獲得的基于ERA5氣壓層產(chǎn)品的ZTD誤差月平均RMS，月平均RMS計(jì)算方法參考式(4)，其中N=720，為月樣本點(diǎn)個(gè)數(shù)(每小時(shí)一個(gè)樣本點(diǎn)，每月按30 d計(jì))。完整的樣本結(jié)構(gòu)如下：

(1) 輸入特征集(氣象特征)包括：地表月平均溫度、地表月平均溫度std、地表月平均水汽、地表月平均水汽std、月平均ZTD、月平均ZTD std、月份。

(2) 輸入特征集(地形特征)包括：230 km高程標(biāo)準(zhǔn)差、4700 km高程標(biāo)準(zhǔn)差、230 km地形切割深度、4700 km地形切割深度、230 km地形起伏度、4700 km地形起伏度、230 km地表粗糙度、4700 km地表粗糙度、230 km高程變異系數(shù)、4700 km高程變異系數(shù)、高程、經(jīng)度。

1.5 數(shù)據(jù)預(yù)處理

(5)

2 模型構(gòu)建

本文以RBF神經(jīng)為基礎(chǔ)構(gòu)建模型。RBF神經(jīng)包含輸入層、隱藏層和輸出層3層。從輸入層到隱含層使用具備徑對(duì)稱特性的激活函數(shù)，實(shí)現(xiàn)非線性變換，從隱含層到輸出層使用線性函數(shù)，實(shí)現(xiàn)線性變換[20]。相較于其他神經(jīng)網(wǎng)絡(luò)，RBF神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)簡(jiǎn)單，訓(xùn)練速度快，逼近能力強(qiáng)，易于解決高維問題[21-22]，但其性能受網(wǎng)絡(luò)結(jié)構(gòu)和規(guī)模影響較大，需合理設(shè)置網(wǎng)絡(luò)參數(shù)和神經(jīng)元個(gè)數(shù)，因此針對(duì)本文所提問題，在2.1節(jié)對(duì)RBF網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行擴(kuò)展，而后針對(duì)網(wǎng)絡(luò)結(jié)構(gòu)規(guī)模問題在2.2節(jié)利用聚類算法確定神經(jīng)元個(gè)數(shù)和激活函數(shù)中心向量，最后2.3節(jié)利用粒子群算法優(yōu)化網(wǎng)絡(luò)權(quán)值。

2.1 擴(kuò)展RBF神經(jīng)網(wǎng)絡(luò)

擴(kuò)展后的RBF神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)如圖3所示，模型由氣象RBF子網(wǎng)絡(luò)、地形RBF子網(wǎng)絡(luò)和線性變換子網(wǎng)絡(luò)3部分組成。氣象和地形RBF子網(wǎng)絡(luò)分別由n和m個(gè)輸入節(jié)點(diǎn)與f和g個(gè)高斯神經(jīng)元組成，線性變換子網(wǎng)絡(luò)由f·g個(gè)線性神經(jīng)元和1個(gè)輸出節(jié)點(diǎn)組成。輸入節(jié)點(diǎn)僅用于傳輸信息，模型運(yùn)算時(shí)輸入樣本被拆分為氣象特征和地形特征兩部分，分別輸入相應(yīng)的RBF子網(wǎng)絡(luò)，RBF子網(wǎng)絡(luò)利用高斯函數(shù)作為基函數(shù)

圖3 擴(kuò)展RBF神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.3 Extended RBF neural network structure

(6)

式中，ui(x)為RBF神經(jīng)元的輸出；x為輸入樣本；ci為高斯函數(shù)中心向量；σi為該神經(jīng)元的標(biāo)準(zhǔn)化常數(shù)。輸入樣本經(jīng)高斯函數(shù)非線性變換后輸出結(jié)果到線性變換子網(wǎng)絡(luò)。線性變換子網(wǎng)絡(luò)中，每一個(gè)氣象RBF神經(jīng)元與每一個(gè)地形RBF神經(jīng)元的輸出各通過一個(gè)線性神經(jīng)元經(jīng)線性加權(quán)輸出結(jié)果到線性變換子網(wǎng)絡(luò)輸出層，輸出節(jié)點(diǎn)對(duì)輸入數(shù)據(jù)線性疊加并輸出結(jié)果。

基于高斯函數(shù)的RBF神經(jīng)元將輸入樣本與高斯函數(shù)中心向量的歐氏距離映射到(0,1)區(qū)間，距離越小，輸出越接近1，反之則接近0，因此模型工作原理是利用RBF神經(jīng)網(wǎng)絡(luò)分別針對(duì)氣象特征和地形特征進(jìn)行分類，將不同大小的ZTD精度認(rèn)定為不同類型的氣象條件和地形環(huán)境共同作用的結(jié)果。

2.2 結(jié)構(gòu)參數(shù)優(yōu)化

氣象RBF子網(wǎng)絡(luò)和地形RBF子網(wǎng)絡(luò)中神經(jīng)元個(gè)數(shù)決定了模型的結(jié)構(gòu)規(guī)模，因而確定合理的神經(jīng)元個(gè)數(shù)對(duì)模型至關(guān)重要，同時(shí)這也是應(yīng)用RBF神經(jīng)網(wǎng)絡(luò)的難點(diǎn)之一[23-24]。本文采用層次聚類法和聚類有效性指標(biāo)分別對(duì)氣象特征集和地形特征集作聚類分析，將最佳聚類數(shù)作為神經(jīng)元個(gè)數(shù)。層次聚類采用歐氏平均距離度量樣本間相似度，采用自下而上的聚合算法生成層次聚類，再通過聚類有效性指標(biāo)選擇最佳聚類數(shù)。聚類有效性指標(biāo)考慮簇內(nèi)數(shù)據(jù)點(diǎn)分布的緊湊度及簇間的分離度，采用式(7)計(jì)算[25]

(7)

由圖5可知，900 ℃下煅燒的生石灰為原料合成的硬硅鈣石明顯可見有絮狀雜質(zhì), 纖維平均直徑約為77 nm。1 000 ℃下煅燒的生石灰為原料合成的硬硅鈣石纖維平均直徑約為82 nm，纖維間搭接規(guī)則，相互交織緊密。1 100 ℃下煅燒的生石灰為原料合成的硬硅鈣石纖維開始出現(xiàn)板結(jié)現(xiàn)象，平均直徑增大到160 nm左右，排布混亂。1 200 ℃下煅燒的生石灰為原料合成的硬硅鈣石纖維板結(jié)現(xiàn)象更加顯著，平均直徑可達(dá)271 nm。

層次聚類還可獲得各聚類中心點(diǎn)，即可作為高斯函數(shù)中心向量，但層次聚類屬于硬聚類，本文所用樣本集事實(shí)上并無嚴(yán)格的類別區(qū)分，不宜使用硬聚類劃分確定聚類中心點(diǎn)。因此，在層次聚類獲得最佳聚類數(shù)的基礎(chǔ)上，利用標(biāo)準(zhǔn)模糊C均值聚類法獲取聚類中心作為高斯函數(shù)中心向量。

2.3 權(quán)值優(yōu)化

模型訓(xùn)練中需調(diào)節(jié)的權(quán)值包括RBF子網(wǎng)絡(luò)高斯函數(shù)標(biāo)準(zhǔn)化常數(shù)和輸出到線性變換子網(wǎng)絡(luò)的連接權(quán)值。本文利用標(biāo)準(zhǔn)粒子群算法搜索最佳權(quán)值，粒子群中每一個(gè)粒子都是權(quán)值空間的一個(gè)解，設(shè)氣象和地形RBF子網(wǎng)絡(luò)神經(jīng)元個(gè)數(shù)分別為f和g，第i個(gè)粒子可表示為Hi=(hi1,hi2,…h(huán)i(f+g+f·g))，粒子i的最佳權(quán)值表示為Pi=(pi1,pi2,…pi(f+g+f·g))，粒子群的最佳權(quán)值表示為Pbest=(bi1,bi2,…bi(f+g+f·g))，粒子群通過式(8)更新速度和權(quán)值

(8)

式中，ω為慣性系數(shù)；s1、s2分別為自學(xué)習(xí)率和群學(xué)習(xí)率；rand1、rand2為[0,1]的隨機(jī)數(shù)。粒子群適應(yīng)度表示為輸出值相對(duì)于目標(biāo)值的均方根。

2.4 數(shù)據(jù)集劃分

由上述模型構(gòu)建過程可知，模型結(jié)構(gòu)規(guī)模、神經(jīng)元數(shù)量和高斯激活函數(shù)中心向量等超參數(shù)分別通過對(duì)訓(xùn)練集進(jìn)行SVD主成分分析降維、層次聚類和標(biāo)準(zhǔn)模糊C均值聚類確定，不必專門配置驗(yàn)證集調(diào)整超參數(shù)，因此將數(shù)據(jù)集劃分為訓(xùn)練集和測(cè)試集兩部分。其中，選擇185個(gè)測(cè)站2017—2019年數(shù)據(jù)樣本作為訓(xùn)練集，185個(gè)測(cè)站2020年數(shù)據(jù)樣本作為測(cè)試集1，余下61個(gè)測(cè)站2017—2020年的數(shù)據(jù)樣本作為測(cè)試集2。測(cè)試集1中185個(gè)測(cè)站地形數(shù)據(jù)參與了模型訓(xùn)練，用于單方面評(píng)估模型通過氣象數(shù)據(jù)估計(jì)ZTD精度能力；測(cè)試集2地形數(shù)據(jù)未參與模型訓(xùn)練，用于完整評(píng)估模型ZTD精度估計(jì)能力。原始數(shù)據(jù)經(jīng)預(yù)處理完成后各數(shù)據(jù)集大小分別為：訓(xùn)練集樣本5821個(gè)、測(cè)試集1樣本2143個(gè)、測(cè)試集2樣本2441個(gè)。

2.5 算法流程

算法流程如圖4所示，具體步驟如下。

圖4 算法流程Fig.4 Algorithm flowchart

(1) 按照1.5節(jié)方法處理樣本數(shù)據(jù)，并劃分?jǐn)?shù)據(jù)集。

(2) 用2.2節(jié)層次聚類和模糊C均值聚類方法處理訓(xùn)練集，獲取神經(jīng)元個(gè)數(shù)和高斯函數(shù)中心值。

(3) 粒子群初始化，確定粒子群規(guī)模、最大訓(xùn)練次數(shù)和相關(guān)常量設(shè)置。

(4) 依據(jù)粒子群權(quán)值和模型結(jié)構(gòu)參數(shù)構(gòu)建并訓(xùn)練神經(jīng)網(wǎng)絡(luò)。

(5) 利用訓(xùn)練結(jié)果結(jié)合訓(xùn)練集的輸出目標(biāo)值求解適應(yīng)度值。

(6) 判斷是否達(dá)到最大訓(xùn)練次數(shù)，結(jié)束訓(xùn)練或依據(jù)適應(yīng)度值更新粒子和群權(quán)值信息并轉(zhuǎn)到步驟(3)。

3 模型性能分析

3.1 模型估計(jì)結(jié)果

訓(xùn)練集層次聚類后獲得氣象特征最佳聚類數(shù)為136，地形特征最佳聚類數(shù)為18，粒子群大小設(shè)置為200，經(jīng)1000次迭代后訓(xùn)練誤差RMS下降到3.7 mm，測(cè)試集1誤差RMS為3.8 mm，測(cè)試集2誤差RMS為4.0 mm。測(cè)試集輸出結(jié)果誤差分布如圖5(a)所示，誤差結(jié)果總體呈正態(tài)分布，測(cè)試集1中82.5%樣本偏差小于5 mm，測(cè)試集2中79.2%樣本偏差小于5 mm。各測(cè)站模型估計(jì)值與參考值對(duì)比如圖5(b)所示，測(cè)試集1中83.8%測(cè)站誤差RMS小于5 mm，測(cè)試集2中82.0%測(cè)站誤差RMS小于5 mm。

圖5 測(cè)試集ZTD精度估計(jì)誤差Fig.5 The ZTD accuracy error for the test set

3.2 結(jié)果分析

測(cè)試集1與訓(xùn)練集中的地形數(shù)據(jù)相同，因此測(cè)試集1估計(jì)結(jié)果反映了模型通過氣象特征估計(jì)ZTD精度的能力，而測(cè)試集2中地形數(shù)據(jù)未參加模型訓(xùn)練，因此測(cè)試集2估計(jì)結(jié)果完整反映了模型在事先未知地形處基于氣象特征和地形特征估計(jì)ZTD精度能力。圖6分別給出了測(cè)試集1和測(cè)試集2中誤差RMS最大和最小的4個(gè)測(cè)站。其中，測(cè)試集1的HLMH和LNYK、測(cè)試集2中的JLYJ和SCGZ誤差RMS較小分別為0.9、1.0、1.6和1.8 mm。由圖6(a)HLMH和LNYK估計(jì)值和參考值對(duì)比可知，在不考慮地形特征情況下模型通過氣象數(shù)據(jù)較好地捕獲了ZTD精度的數(shù)值特征和變化趨勢(shì)，達(dá)到了預(yù)期的估計(jì)效果。由圖6(b)中JLYJ和SCGZ參考值與估計(jì)值對(duì)比可知，模型輸出結(jié)果準(zhǔn)確估計(jì)了ZTD精度大小和周期變化趨勢(shì)(冬季精度較高，夏季精度較低，變化周期為12個(gè)月)，模型不僅對(duì)已訓(xùn)練過測(cè)站具備較好的估計(jì)能力，而且對(duì)事先未接觸過的測(cè)站也表現(xiàn)出良好的估計(jì)能力，說明模型從數(shù)據(jù)集中成功提取了相應(yīng)特征，具備較好的泛化能力。

圖6 測(cè)試集1和測(cè)試集2中ZTD精度誤差RMS最大和最小的4個(gè)GNSS測(cè)站Fig.6 The four GNSS stations with maximum and minimum ZTD accuracy error RMS in test set 1 and test set 2

但同時(shí)也注意到，測(cè)試集中還存在少數(shù)誤差RMS較大的測(cè)站，如圖6中測(cè)試集1的XJSH和XJDS(11.1 mm、12.7 mm)和測(cè)試集2的YNMH和XZBG(7.2 mm、6.4 mm)。由圖6可知，模型基本捕捉到了上述測(cè)站ZTD精度變化的周期性趨勢(shì)，但對(duì)ZTD精度的數(shù)值特征估計(jì)偏差較大，測(cè)站XJDS和XJSH估計(jì)值整體偏小，而測(cè)站XZBG的模型估計(jì)值則是先偏小后偏大，原因有多方面：一是訓(xùn)練集中地形數(shù)據(jù)受GNSS測(cè)站數(shù)量限制，特征提取和分類還不夠完備；二是參考值數(shù)據(jù)質(zhì)量問題，現(xiàn)有參考文獻(xiàn)中對(duì)陸態(tài)網(wǎng)部分測(cè)站對(duì)比分析結(jié)果表明，其數(shù)據(jù)精度在4 mm左右[14]，但現(xiàn)有比對(duì)分析還不夠全面，可能存在部分測(cè)站精度較差的情況，如測(cè)站XJDS，其ZTD精度參考值發(fā)生階梯狀變化很可能就是測(cè)站自身硬件或軟件升級(jí)改善引起的。

進(jìn)一步對(duì)比觀察圖6中不同測(cè)站估計(jì)值和參考值可以發(fā)現(xiàn)，估計(jì)結(jié)果誤差RMS與ZTD精度估計(jì)值大小存在一定關(guān)系，誤差RMS較大的測(cè)站ZTD精度估計(jì)值也較小，據(jù)此分析了輸出結(jié)果誤差絕對(duì)值與輸出結(jié)果的關(guān)系，如圖7所示，圖中每個(gè)點(diǎn)代表一個(gè)測(cè)試集樣本，ZTD精度估計(jì)結(jié)果小于10 mm時(shí)，誤差整體較小，而當(dāng)樣本估計(jì)結(jié)果大于10 mm時(shí)，出現(xiàn)少數(shù)樣本誤差較大的情況。因此在使用本模型估計(jì)ZTD精度時(shí)，當(dāng)估計(jì)值小于10 mm時(shí)具有較高可信度，而大于10 mm時(shí)則需注意存在偏差較大的可能性。

考慮到作為參考值的陸態(tài)網(wǎng)ZTD產(chǎn)品本身精度為4 mm左右，本文認(rèn)為估值誤差小于5 mm的結(jié)果是合理有效的，結(jié)合圖6和圖7，可以認(rèn)為在絕大部分地區(qū)，模型估計(jì)結(jié)果是有效的。

圖7 測(cè)試集輸出結(jié)果與誤差絕對(duì)值分布Fig.7 The model output results and absolute error distribution of test set

采用本文構(gòu)建模型，以2019年3月和8月ERA5氣象數(shù)據(jù)為基礎(chǔ)，按照1°×1°水平間隔計(jì)算了中國(guó)內(nèi)陸及其周邊區(qū)域ZTD精度，結(jié)果如圖8所示，東北及華北平原ZTD精度較高，青藏高原和云貴高原次之，四川盆地和云貴交接處ZTD精度偏低，新疆塔里木盆地和河西走廊及黃土高原沿線以北等地區(qū)精度最低；長(zhǎng)江中下游以南地區(qū)精度大小與季節(jié)密切相關(guān)，3月較小，8月較大。由上述結(jié)果可知，基于NWM模型獲取的ZTD在地勢(shì)平坦地區(qū)精度較高，而地形變化劇烈區(qū)域精度較低；內(nèi)陸地區(qū)季節(jié)性變化較弱，沿海地區(qū)季節(jié)性變化明顯。另須特別注意，本文數(shù)據(jù)集中缺少海上測(cè)站相關(guān)信息，因此基于本文模型生成的海上ZTD精度不具備參考價(jià)值。

圖8 2019年3月和8月的ZTD精度分布Fig.8 The ZTD accuracy distribution for March and August in 2019

4 結(jié) 論

本文以NWM氣象數(shù)據(jù)和地形特征數(shù)據(jù)為特征向量，采用層次聚類法和模糊C均值聚類確定網(wǎng)絡(luò)規(guī)模，使用粒子群算法優(yōu)化模型參數(shù)，構(gòu)建了耦合粒子群算法與擴(kuò)展RBF神經(jīng)網(wǎng)絡(luò)的ZTD精度估計(jì)模型。為驗(yàn)證模型有效性，以ERA5為NWM特例，進(jìn)行了模型訓(xùn)練和結(jié)果驗(yàn)證。研究結(jié)果表明：模型可在絕大部分區(qū)域有效估計(jì)ZTD月平均精度，實(shí)現(xiàn)不依賴于外部參考基準(zhǔn)的ZTD精度估計(jì)，可在未配置GNSS測(cè)站的區(qū)域提供有參考價(jià)值的估值結(jié)果。同時(shí)也看到，本文所用數(shù)據(jù)集規(guī)模有限，測(cè)站數(shù)量較小，下一步將研究把數(shù)據(jù)集規(guī)模擴(kuò)展到全球，進(jìn)一步提升并驗(yàn)證模型性能。