云貴川地區(qū)ZTD垂直剖面格網(wǎng)模型

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

1 桂林理工大學(xué)測繪地理信息學(xué)院，桂林市雁山街319號(hào)，541006 2 湖北科技學(xué)院資源環(huán)境科學(xué)與工程學(xué)院，湖北省咸寧市咸寧大道88號(hào)，437100

由于對流層天頂延遲(ZTD)在天頂方向上的誤差約為2 m，隨高度角的降低誤差會(huì)增大至20 m[1]，因此對流層天頂延遲改正對全球?qū)Ш叫l(wèi)星系統(tǒng)GNSS實(shí)時(shí)高精度導(dǎo)航定位及GNSS水汽反演具有重要意義。目前，對流層天頂延遲模型主要分為2大類：1)需要實(shí)測氣象參數(shù)的模型，如Hopfield模型[2]和Saastamoinen模型[3]；2)實(shí)時(shí)模型，如UNB3m模型[4]和EGNOS模型[5]。但上述2類模型均存在一些缺陷：用戶位置的實(shí)時(shí)氣象參數(shù)難以獲取，限制了氣象參數(shù)模型的應(yīng)用；UNB3m模型和EGNOS模型的時(shí)空分辨率較低，限制了實(shí)時(shí)模型的應(yīng)用。GPT2w與GPT3模型是全球格網(wǎng)經(jīng)驗(yàn)?zāi)Ｐ停珿PT3模型改善了GPT2w模型中的映射函數(shù)體系，性能優(yōu)異且應(yīng)用廣泛，但不能提供對流層天頂延遲的垂直剖面信息。近年來，眾多學(xué)者利用大氣再分析資料對對流層垂直剖面函數(shù)展開廣泛研究[6]，采用二次多項(xiàng)式[7]和負(fù)指數(shù)函數(shù)[8-10]進(jìn)行模型構(gòu)建。此外，姚宜斌等[11]利用多源數(shù)據(jù)建立考慮季節(jié)變化的對流層天頂延遲模型；毛健等[12]、黃良珂等[13]建立同時(shí)顧及高程和季節(jié)變化的對流層天頂延遲模型。上述模型均表現(xiàn)出較好的性能。

在地形起伏極大、氣候復(fù)雜多變、晝夜大氣差異顯著的云貴川地區(qū)，已有模型難以提供該地區(qū)實(shí)時(shí)高精度對流層天頂延遲信息。基于此，本文利用2015～2017年高時(shí)空分辨率的ERA5再分析資料積分計(jì)算分層ZTD，分析高程尺度因子的精細(xì)時(shí)間變化特征，采用負(fù)指數(shù)函數(shù)形式構(gòu)建顧及高程尺度因子精細(xì)時(shí)空變化的云貴川地區(qū)ZTD垂直剖面格網(wǎng)模型(YZTD-H模型)。其中，YZTD-H1模型顧及了年周期、半年周期、日周期；YZTD-H2模型在YZTD-H1模型的基礎(chǔ)上增加了半日周期。本文還利用2018年云貴川地區(qū)探空站分層ZTD數(shù)據(jù)對模型進(jìn)行精度檢驗(yàn)。

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

ERA5大氣再分析資料(https:∥apps.ecmwf.int/datasets/data/interim-full-daily)是歐洲中尺度天氣預(yù)報(bào)中心ECMWF第五代全球氣候再分析數(shù)據(jù)集，可提供逐小時(shí)的地表參數(shù)和等壓分層數(shù)據(jù)，有望在未來得到廣泛應(yīng)用。ERA5數(shù)據(jù)的水平分辨率為0.25°×0.25°，垂直分辨率為37層。探空站數(shù)據(jù)來源于美國俄懷明大學(xué)官網(wǎng)(http:∥weather.uwyo.edu/upperair/sounding.html)，該數(shù)據(jù)可提供從地面到近地空間約30 km的氣壓、溫度、風(fēng)速、相對濕度等氣象參數(shù)的分層數(shù)據(jù)和大氣水汽等地表參數(shù)，每日進(jìn)行2次探測(UTC00:00和12:00)。探空站數(shù)據(jù)為實(shí)測值，具有較高的精度和可靠性，可作為模型檢驗(yàn)的參考值。本文選取云貴川地區(qū)9個(gè)探空站對YZTD-H模型進(jìn)行檢驗(yàn)。

2 模型建立

覃澤穎等[14]對ERA5數(shù)據(jù)計(jì)算得到的ZTD精度進(jìn)行評估，結(jié)果表明，ERA5數(shù)據(jù)可作為ZTD模型建立的數(shù)據(jù)源。本文利用2015～2017年云貴川地區(qū)的ERA5數(shù)據(jù)，通過積分法計(jì)算得到格網(wǎng)垂直剖面ZTD。

2.1 高程尺度因子的特性分析

負(fù)指數(shù)函數(shù)在表示ZTD垂直剖面變化方面具有較好的優(yōu)越性，公式如下：

(1)

式中，ZTDt為目標(biāo)高程Ht處的ZTD，ZTDr為起算高程Hr處的ZTD，Hz為ZTD的高程尺度因子。利用2015～2017年云貴川地區(qū)ERA5格網(wǎng)垂直剖面ZTD數(shù)據(jù)對Hz進(jìn)行擬合，得到每個(gè)格網(wǎng)點(diǎn)的逐小時(shí)Hz。為分析Hz的精細(xì)時(shí)間變化特性，選取云貴川地區(qū)的2個(gè)格網(wǎng)點(diǎn)，計(jì)算這2個(gè)格網(wǎng)點(diǎn)2015～2017年的日均Hz，對其進(jìn)行年周期和半年周期的三角函數(shù)擬合，并采用快速傅里葉變換法分析其周期性，結(jié)果如圖1所示。圖1(a)、(b)分別為2個(gè)格網(wǎng)點(diǎn)3 a日均Hz的時(shí)間序列及其擬合線；圖1(c)、(d)分別為2個(gè)格網(wǎng)點(diǎn)Hz的頻譜分析結(jié)果。

圖1 Hz的季節(jié)變化及其頻譜分析Fig.1 Seasonal variation of Hz and its spectral analysis

由圖1可見，Hz存在明顯的年周期和半年周期，其中格網(wǎng)點(diǎn)1的半年周期較格網(wǎng)點(diǎn)2更加明顯。為分析Hz的精細(xì)時(shí)間變化特性，同樣選取上述2個(gè)格網(wǎng)點(diǎn)計(jì)算年均UTC 00:00～23:00的逐小時(shí)Hz，對其進(jìn)行日周期和半日周期的三角函數(shù)擬合，并采用快速傅里葉變換法分析其周期性(圖2)。由圖可見，Hz存在顯著的日周期和半日周期。

圖2 Hz的日變化及其頻譜分析Fig.2 Daily variation of Hz and its spectral analysis

2.2 模型建立

在上述Hz周期特性分析的基礎(chǔ)上，對Hz進(jìn)行顧及年均值、年周期、半年周期、日周期、半日周期的三角函數(shù)擬合，表達(dá)式如下：

(2)

由圖3可見，Hz的各擬合系數(shù)存在顯著的區(qū)域性。由于云貴川地區(qū)的地勢特征為西北高、東南低，因此各擬合系數(shù)呈現(xiàn)出由東南向西北逐漸減小的趨勢。多個(gè)擬合系數(shù)在四川盆地和西北高海拔地區(qū)存在顯著差異，說明擬合系數(shù)受海拔影響較大。擬合系數(shù)在東南沿海與內(nèi)陸地區(qū)也存在顯著差異，可能是沿海與內(nèi)陸的氣候差異和晝夜大氣差異所致。由此可知，Hz具有精細(xì)的時(shí)空分布特征，需要對每個(gè)ERA5格網(wǎng)點(diǎn)的Hz進(jìn)行精細(xì)的周期性擬合。考慮擬合系數(shù)存儲(chǔ)量及半日周期對Hz的影響程度后，本文選用2種方式對2015～2017年云貴川地區(qū)ERA5格網(wǎng)Hz進(jìn)行擬合：顧及年周期、半年周期、日周期的擬合以及顧及年周期、半年周期、日周期、半日周期的擬合。Hz的模型系數(shù)以0.25°×0.25°水平分辨率的格網(wǎng)進(jìn)行存儲(chǔ)，最終構(gòu)建出云貴川地區(qū)ZTD垂直剖面格網(wǎng)模型YZTD-H1和YZTD-H2。用戶只需提供經(jīng)緯度、高程、年積日及UTC時(shí)等數(shù)據(jù)，模型便可結(jié)合式(1)和式(2)將用戶位置從起算高程處的ZTD改正到目標(biāo)高程處。

圖3 Hz各擬合系數(shù)的分布Fig.3 Distribution of fitting coefficients of Hz

3 模型精度驗(yàn)證

3.1 模型精度對比

利用云貴川地區(qū)9個(gè)探空站數(shù)據(jù)對YZTD-H1和YZTD-H2模型進(jìn)行精度檢驗(yàn)，并將其與性能優(yōu)異的GPT2w、GPT3模型進(jìn)行精度對比(GPT2w模型和GPT3模型的水平分辨率均為1°×1°)。UNB3模型可將天頂靜力學(xué)延遲(zenith hydrostatic delay, ZHD)和天頂濕延遲(zenith wet delay, ZWD)由海平面高度改正至目標(biāo)高度，但該模型需要?dú)庀髤?shù)，而GPT2w與GPT3模型可提供全球任意位置的氣象參數(shù)，因此將二者結(jié)合即可實(shí)現(xiàn)ZTD的垂直剖面改正。由于作為參考值的探空站數(shù)據(jù)并非從海平面開始分層，因此需要對UNB3模型進(jìn)行推導(dǎo)，使其能夠?qū)θ我飧叨鹊腪HD和ZWD進(jìn)行高程改正。推導(dǎo)結(jié)果如式(3)、式(4)所示：

(3)

(4)

式中，ZHDt和ZHDr分別為目標(biāo)高度和起算高度的ZHD值；ZWDt和ZWDr分別為目標(biāo)高度和起算高度的ZWD值；Ht和Hr分別為目標(biāo)高度和起算高度；β為溫度梯度；T0為熱力學(xué)溫度；Rd為干氣體常量(Rd=287.053 8 J/kg)；g為地表重力加速度；λ′=1+λ，λ為水汽梯度。首先采用積分法計(jì)算2018年云貴川地區(qū)探空站分層ZTD值及首層ZHD/ZWD值，將分層ZTD值作為模型檢驗(yàn)的參考值；然后以探空站首層高度作為起算高度，以探空站首層ZTD作為起算值、其他層高作為目標(biāo)高度，利用YZTD-H1、YZTD-H2模型將探空站首層ZTD依次改正至其他各層。同樣，首先以探空站首層高度作為起算高度，以探空站首層ZHD/ZWD作為起算值、其他層高作為目標(biāo)高度，利用GPT2w、GPT3模型將探空站首層ZHD/ZWD依次改正至其他各層，然后對目標(biāo)高度的ZHD和ZWD求和，即可得到目標(biāo)高度的ZTD值。將偏差bias和均方根誤差RMSE作為精度指標(biāo)，計(jì)算2018年云貴川地區(qū)9個(gè)探空站UTC 00:00和12:00的各模型分層ZTD值，得到各模型的總體bias和RMSE，結(jié)果如表1所示(單位mm)。

由表1可見，YZTD-H1和YZTD-H2模型均表現(xiàn)出接近于0的平均bias值，分別為-7.81 mm和-7.60 mm，其中YZTD-H2模型略優(yōu)于YZTD-H1模型；GPT2w和GPT3模型均表現(xiàn)出較大的正bias值，說明二者存在較大的系統(tǒng)偏差，其中GPT3模型略優(yōu)于GPT2w模型。YZTD-H1和YZTD-H2模型的平均RMSE分別為50.50 mm和50.44 mm，YZTD-H2模型略優(yōu)于YZTD-H1模型；GPT2w和GPT3模型的平均RMSE分別為79.94 mm和79.93 mm，GPT3模型略優(yōu)于GPT2w模型。由此可知，YZTD-H1和YZTD-H2模型優(yōu)于GPT2w和GPT3模型：YZTD-H1模型較GPT2w和GPT3模型的RMSE分別減少29.44 mm和29.43 mm；YZTD-H2模型較GPT2w和GPT3模型的RMSE分別減少29.50 mm和29.49 mm。綜上可知，YZTD-H1和YZTD-H2模型具有更高的精度和更好的穩(wěn)定性。

表1 不同模型的精度統(tǒng)計(jì)

3.2 各模型在探空站的精度對比

為分析各模型在探空站的精度表現(xiàn)，計(jì)算4種模型在各探空站的年均bias和RMSE，結(jié)果如圖4所示。

圖4 不同模型的精度分布Fig.4 Precision distribution of different models

由圖4可見，各模型的年均bias和RMSE存在一定的區(qū)域性，可能與云貴川地區(qū)復(fù)雜的地形、多變的氣候和較大的晝夜大氣差異有關(guān)。YZTD-H1和YZTD-H2模型在各探空站的年均bias總體接近于0，緯度較低的2個(gè)探空站表現(xiàn)出較大的負(fù)偏差；GPT2w和GPT3模型在各探空站的年均bias總體呈現(xiàn)出較大的正偏差。YZTD-H1和YZTD-H2模型在各探空站的年均RMSE約為50 mm，東北方向的探空站精度更高；GPT2w和GPT3模型在各探空站的年均RMSE均大于60 mm，西南方向的探空站精度更高。由此可知，YZTD-H1和YZTD-H2模型對地理位置的適應(yīng)性更強(qiáng)。

3.3 模型精度的季節(jié)變化特征

由于顧及半日周期的YZTD-H2模型與未顧及半日周期的YZTD-H1模型性能相當(dāng)，且YZTD-H2模型的RMSE比YZTD-H1模型僅高了0.06 mm，因此選用YZTD-H1模型能夠在保證模型精度的前提下減少2個(gè)模型系數(shù)，減小模型計(jì)算量。由于GPT2w模型與GPT3模型精度相當(dāng)，因此本文僅對YZTD-H1模型與GPT3模型進(jìn)行精度對比。為分析模型精度的季節(jié)變化特征，利用YZTD-H1模型與GPT3模型計(jì)算2018年云貴川地區(qū)各探空站的日均bias和RMSE，結(jié)果如圖5所示。

圖5 不同模型日均bias和RMSE的散點(diǎn)分布Fig.5 Scatter distribution of daily bias and RMSE for different models

由圖5可見，GPT3模型的日均bias呈現(xiàn)出明顯的季節(jié)性，夏季的正偏差較大；YZTD-H1模型的日均bias均在略小于0的附近波動(dòng)。GPT3模型的日均RMSE同樣呈現(xiàn)出明顯的季節(jié)性，夏季的RMSE大于春、冬季節(jié)，可能是因?yàn)橄募咎厥獾臍夂驅(qū)е麓髿馑推渌麣庀髤?shù)發(fā)生變化，進(jìn)而導(dǎo)致模型難以精準(zhǔn)預(yù)測ZTD；YZTD-H1模型的日均RMSE穩(wěn)定在50 mm附近。由此可知，YZTD-H1模型能夠在氣候條件復(fù)雜和晝夜大氣差異大的云貴川地區(qū)表現(xiàn)出更好的適應(yīng)性。

3.4 模型在不同插值高度上的精度變化

為分析模型在不同插值高度上的精度變化，將探空站首層高度與其他層的高度差分為6個(gè)區(qū)間，統(tǒng)計(jì)YZTD-H1模型和GPT3模型在各區(qū)間內(nèi)的bias和RMSE，結(jié)果如圖6所示。

圖6 不同模型bias和RMSE隨插值高度的變化Fig.6 Variations with interpolated height of bias and RMSE for different models

由圖6可見，GPT3模型的bias在各插值高度上均表現(xiàn)出正偏差，且bias隨插值高度的增大而減小；YZTD-H1模型的bias在4 km以下為負(fù)偏差，在4 km以上為正偏差，整體偏差較小。GPT3模型的RMSE隨插值高度的增大而減小；YZTD-H1模型的RMSE隨插值高度的增大呈現(xiàn)出先增大后減小的變化趨勢，整體變化較為平緩。8 km以上的GPT3模型精度與YZTD-H1模型相當(dāng)，說明GPT3模型只有在較高的插值高度上才具有較好的適應(yīng)性，而YZTD-H1在任意插值高度上均有很好的適應(yīng)性。

3.5 模型在不同誤差區(qū)間內(nèi)的分布特征

為分析模型在不同誤差區(qū)間內(nèi)的分布特征，將bias和RMSE均分為7個(gè)區(qū)間，分別統(tǒng)計(jì)YZTD-H1模型和GPT3模型的日均bias和RMSE在各區(qū)間的數(shù)量，結(jié)果如圖7所示。

圖7 不同模型日均bias和RMSE在不同區(qū)間內(nèi)的分布Fig.7 Distribution of daily bias and RMSE for different models in different interval

由圖7可見，GPT3模型的bias和RMSE在各區(qū)間內(nèi)均勻分布，說明該模型穩(wěn)定性較差；YZTD-H1模型的bias集中分布在0 附近，RMSE集中分布在40～60 mm之間，進(jìn)一步說明YZTD-H1模型具有更優(yōu)的性能與更好的穩(wěn)定性。

4 結(jié) 語

1)YZTD-H1模型和YZTD-H2模型的總體精度(年均RMSE)均優(yōu)于GPT2w模型和GPT3模型，YZTD-H2模型精度略優(yōu)于YZTD-H1模型，二者在各探空站的年均bias和RMSE表現(xiàn)穩(wěn)定，且并未表現(xiàn)出較為明顯的區(qū)域性。

2)在YZTD-H1模型與GPT3模型的不同季節(jié)精度對比中，GPT3模型表現(xiàn)出明顯的季節(jié)性，而YZTD-H1模型的全年性能表現(xiàn)穩(wěn)定。

3)在YZTD-H1模型與GPT3模型的不同插值高度精度對比中，GPT3模型受插值高度的影響較大，而YZTD-H1模型在各插值高度區(qū)間內(nèi)表現(xiàn)穩(wěn)定，能將RMSE控制在50 mm左右，說明負(fù)指數(shù)函數(shù)形式能夠較好地將垂直剖面ZTD模型化。

4)在YZTD-H1模型與GPT3模型的不同誤差區(qū)間對比中，YZTD-H1模型的日均RMSE集中分布在40～60 mm之間，而GPT3模型在各區(qū)間內(nèi)均勻分布，進(jìn)一步說明YZTD-H1模型具有更優(yōu)的性能與更好的穩(wěn)定性。

綜上所述，YZTD-H1模型在不同地理位置、不同季節(jié)、不同插值高度上均表現(xiàn)出良好的穩(wěn)定性和優(yōu)異的性能，更適用于云貴川地區(qū)ZTD的垂直剖面插值計(jì)算。本研究對云貴川地區(qū)GNSS實(shí)時(shí)高精度導(dǎo)航定位和GNSS水汽反演具有很好的參考價(jià)值。