基于Emardson模型的GNSS-PWV轉(zhuǎn)換系數(shù)精化

李黎， 侯曉玲， 張?chǎng)?高穎， 趙偉， 陳國(guó)棟

1 蘇州科技大學(xué)地理科學(xué)與測(cè)繪工程學(xué)院， 蘇州 215009 2 蘇州科技大學(xué)北斗導(dǎo)航與環(huán)境感知研究中心， 蘇州 215009

0 引言

水汽在大氣中所占比例很小，僅占大氣體積的0.1%～0.3%，但它卻是龍卷風(fēng)、雷暴雨和臺(tái)風(fēng)等災(zāi)害性天氣形成的重要因素之一(宋淑麗等, 2004; 畢研盟等, 2006)，獲取高時(shí)空分辨率的水汽資料可以提高降水預(yù)報(bào)的精度，減少災(zāi)害性天氣給國(guó)民帶來的經(jīng)濟(jì)損失.由于傳統(tǒng)水汽觀測(cè)方法的時(shí)空分辨率較低，限制了降水預(yù)報(bào)精度，實(shí)踐中使用并不廣泛(Baker et al., 2001).自Bevis等(1992)提出GNSS氣象學(xué)以后，GNSS技術(shù)以其全球覆蓋、全天候、精度高、垂直分辨率高(≤1 km)和費(fèi)用低等優(yōu)點(diǎn)逐漸成為獲取大氣可降水量 (Precipitable Water Vapour, PWV) 的重要手段之一.由GNSS技術(shù)(李黎等, 2012)解算得到的對(duì)流層延遲(Zenith Total Delay, ZTD)減去根據(jù)Saastamoinen模型計(jì)算的對(duì)流層干延遲(Zenith Hydrostatic Delay, ZHD)可得對(duì)流層濕延遲(Zenith Wet Delay, ZWD)，ZWD再乘以水汽轉(zhuǎn)換系數(shù)(K值)可得到GNSS-PWV，可見K值是影響GNSS-PWV反演精度的主要因素之一，因此研究并建立高精度的K值模型非常必要.

近40年來，GNSS氣象學(xué)的發(fā)展開拓了高時(shí)空分辨率GNSS-PWV的應(yīng)用研究(Bevis et al., 1994; Wen et al., 2012).而K值的精度在一定程度上受限于大氣加權(quán)平均溫度(Tm)的精度，不少學(xué)者已深入研究Tm并建立區(qū)域Tm模型.計(jì)算Tm需用到地表溫度、氣壓和高程等信息，這些氣象要素的實(shí)測(cè)精度和實(shí)時(shí)程度會(huì)進(jìn)一步影響Tm的最終精度.Emardson和Derks(2000)基于歐洲數(shù)十個(gè)探空站資料，建立一個(gè)與位置和年積日相關(guān)的K值計(jì)算模型(Emardson模型).該模型未引入Tm參數(shù)，從而避免氣象要素的誤差積累，進(jìn)一步簡(jiǎn)化了計(jì)算過程，在實(shí)時(shí)反演GNSS-PWV中有著廣泛的應(yīng)用前景.鑒于Emardson模型在不同地區(qū)適用性不同，姚朝龍等(2015)建立了一個(gè)適用于低緯度地區(qū)的K值模型，發(fā)現(xiàn)其精度主要取決于緯度和高程，與單站模型或者通過Tm計(jì)算的K值精度相當(dāng)；劉立龍等(2016)選取新疆地區(qū)9個(gè)無線電探空站連續(xù)5年的探空資料，根據(jù)積分法獲取的K值分析Emardson模型在該地區(qū)的時(shí)空適用性；劉立龍等(2017)分析了地理位置對(duì)于K值的影響，在廣西地區(qū)建立的本地化K值模型顯示顧及高程的Emardson模型精度更好；黃良珂等(2019)利用不同的K值模型反演拉薩站的GNSS-PWV，其結(jié)果與探空PWV具有較好的一致性，表明顧及高程的Emardson-H模型適用于缺乏氣象參數(shù)的高海拔地區(qū)；侯曉玲等(2021)根據(jù)長(zhǎng)三角地區(qū)7個(gè)探空站2年的探空數(shù)據(jù)，建立了顧及與不顧及高程的兩個(gè)精化模型，發(fā)現(xiàn)在低海拔地區(qū)兩種模型得到的GNSS-PWV精度相當(dāng).

鑒于Emardson模型具有地區(qū)差異性，中國(guó)地區(qū)缺乏統(tǒng)一的高精度K值預(yù)報(bào)模型，本文利用2018年中國(guó)地區(qū)的42個(gè)探空站資料，建立顧及與不顧及高程的兩套Emardson精化模型，并以2019年K值作為參考，分析以上兩個(gè)模型的K值預(yù)報(bào)精度，評(píng)估其在中國(guó)地區(qū)的適用性，以期快速獲取GNSS-PWV，提高其在氣象預(yù)報(bào)領(lǐng)域的實(shí)時(shí)應(yīng)用效率.

1 數(shù)據(jù)處理方法與數(shù)據(jù)來源及分析

1.1 K值計(jì)算方法

利用探空數(shù)據(jù)所提供的氣壓、水汽壓和溫度等氣象參數(shù)和積分法是目前應(yīng)用最廣泛的K值計(jì)算方法.

(1)

(2)

式中，Pv為平均水汽壓(hPa)，T為絕對(duì)溫度(K)，Δh為大氣的厚度(m).T，Δh均可在探空資料中得到，Tm是指探空站上空水汽壓與絕對(duì)溫度沿天頂方向的積分值，是計(jì)算K值的關(guān)鍵參數(shù).大氣水汽一般分布在離地面12 km以下，故可以通過數(shù)值積分法來逼近Tm真值.在數(shù)值積分法中，觀測(cè)誤差與逼近誤差的綜合影響小于1 K(劉焱雄等, 2000)，式(2)是目前公認(rèn)精度最高的Tm計(jì)算方法(Duan et al., 1996; Ingold et al., 1998).Pv無法直接觀測(cè)，但能通過露點(diǎn)溫度和飽和水汽壓公式間接計(jì)算,

Pv=6.112exp[17.62td/(243.12+td)],

(3)

式中，td為露點(diǎn)溫度(℃).

1.2 PWV計(jì)算方法

根據(jù)GNSS技術(shù)解算，可以求得GNSS-ZTD(精度優(yōu)于6 mm)，一般時(shí)間分辨率為5 min.根據(jù)地面緯度、氣壓、高程等信息，基于Saastamoinen模型可精確計(jì)算ZHD(Saastamoinen, 1972)：

(4)

式中，Ps為表面壓力(hPa)，φ為緯度，H為大地高度(km).測(cè)站位置確定后，ZHD是地表壓力的函數(shù)，壓力每增加1 hPa只會(huì)引起大約0.2 mm的ZHD誤差(Tregoning and Herring, 2006).

通過GNSS-ZTD減去ZHD，可得到ZWD，ZWD與K相乘可得到GNSS-PWV.

(5)

1.3 數(shù)據(jù)來源及其分析

1.3.1 數(shù)據(jù)來源與測(cè)站信息

本文的探空資料來源于2018—2019年全國(guó)均勻分布的42個(gè)無線電探空站每天2次的實(shí)測(cè)數(shù)據(jù)，采樣間隔為12 h，垂直分辨率為1 km.實(shí)測(cè)數(shù)據(jù)可從美國(guó)懷俄明大學(xué)天氣數(shù)據(jù)網(wǎng)(http:∥weather.uwyo.edu/upperair/sounding.html)獲取，包括氣壓、露點(diǎn)溫度、高度、地面溫度、相對(duì)濕度等氣象要素，該機(jī)構(gòu)以世界氣象組織制定的規(guī)則為依據(jù)，已對(duì)資料進(jìn)行質(zhì)量檢查.GNSS數(shù)據(jù)可從中國(guó)地震局GNSS產(chǎn)品服務(wù)平臺(tái)(http:∥www.cgps.ac.cn/index.htm)獲取(Zhao et al., 2015; Yu et al., 2019)，其中包括GNSS-PWV及其梯度、氣壓、溫度等參數(shù)，時(shí)間分辨率為1 h.

本文將用2018年中國(guó)地區(qū)的42個(gè)探空站積分法所得的K值建立兩套Emardson模型，之后用2019年的42個(gè)探空站和11個(gè)GNSS站分別驗(yàn)證兩套預(yù)報(bào)模型的K值精度及其PWV反演精度.圖1是中國(guó)地區(qū)42個(gè)探空站和11個(gè)GNSS站的分布圖，表1和表2分別是探空站和GNSS站的地理位置信息(包括站名、站號(hào)、代碼、經(jīng)緯度和高程等).

1.3.2K值與測(cè)站位置的相關(guān)性分析

我國(guó)地處太平洋西岸，橫跨多個(gè)時(shí)區(qū)，本文以式(1)計(jì)算的2018年42個(gè)探空站K值作為真值，分析研究K值與測(cè)站經(jīng)度和緯度之間的相關(guān)性，結(jié)果見圖2.由圖2可知，K值與測(cè)站經(jīng)度之間的相關(guān)系數(shù)為0.592，二者的相關(guān)性較小；而K值與測(cè)站緯度之間的相關(guān)系數(shù)為0.735，變化趨勢(shì)呈線性負(fù)相關(guān)，即緯度越高，測(cè)站K值則越小.因此，在建立K值精化模型時(shí)，測(cè)站經(jīng)度的影響可以不予考慮，但需要把緯度作為建模的影響因子之一.

中國(guó)地區(qū)的地形總體呈現(xiàn)西高東低，且起伏較大，因此有必要分析K值與測(cè)站高程之間的線性關(guān)系.由圖3可知，K值與測(cè)站高程呈負(fù)相關(guān)，K值隨著測(cè)站高程的增加而逐漸下降，兩者之間的相關(guān)系數(shù)為0.941，尤其高程在0～1000 m區(qū)間時(shí)，下降趨勢(shì)更為明顯.因此，在中國(guó)地區(qū)建立K值精化模型時(shí)，尤其需要考慮測(cè)站高程對(duì)K值模型精度的影響，要將高程作為建模的主要影響因子之一.

圖1 中國(guó)地區(qū)探空站及GNSS站地理位置分布圖Fig.1 Distribution map of selected radiosondes and GNSS stations in China

表1 中國(guó)地區(qū)42個(gè)探空站地理位置Table 1 Geographical position of 42 selected radiosonde stations in China

表2 中國(guó)地區(qū)11個(gè)GNSS站地理位置Table 2 Geographical position of 11 selected GNSS stations in China

圖2 K值與探空站經(jīng)度和緯度之間的相關(guān)性分析Fig.2 Correlation analysis between K value and longitude and latitude of the radiosonde stations

圖3 K值與探空站高程的相關(guān)性分析Fig.3 Correlation analysis between K value and altitude of radiosonde stations

2 中國(guó)地區(qū)Emardson模型的建立及其精度分析

2.1 Emardson模型的建立

根據(jù)年積日、測(cè)站緯度和高程等時(shí)空特性與K值之間的相關(guān)性，分別建立基于測(cè)站緯度和年積日的Emardson-I模型和基于測(cè)站緯度、高程和年積日的Emardson-H模型(Emardson and Derks, 2000).

(6)

+α4H,

(7)

式中，φ為大地緯度(°)，doy為年積日，H為高程(m)，α0、α1、α2、α3、α4為模型系數(shù)，α0為常數(shù)，α1為φ的系數(shù)，α4為H的系數(shù).

本文將2018年中國(guó)地區(qū)42個(gè)探空數(shù)據(jù)計(jì)算的K值作為真值，使用多元線性擬合法解算得到模型系數(shù)(如表3).

表3 中國(guó)地區(qū)的Emardson-I和Emardson-H模型系數(shù)Table 3 The model coefficients of Emardson-I and Emardson-H in China

結(jié)合表3、式(6)和式(7)可知，α1代表K值與緯度之間的關(guān)系，負(fù)值表示K值隨探空站緯度的升高而減小，緯度每升高1°，K值約減小0.00023；α4代表K值與高程之間的關(guān)系，負(fù)值表示K值隨探空站處高程的增加而減小，高程每增加1 m，K值則減小0.0000012.因此，在地形平緩的低海拔地區(qū)，K值隨高程變化的影響較小；而在地形起伏較大時(shí)，尤其是我國(guó)西部3000 m以上的高海拔地區(qū)，高程對(duì)K值的影響不可忽視.

2.2 K值預(yù)報(bào)精度分析

利用2.1節(jié)得到的Emardson-I和Emardson-H模型預(yù)報(bào)2019年42個(gè)探空站處的K值，以2019年積分法K值作為參考真值，評(píng)估兩套K值模型的預(yù)報(bào)精度.限于篇幅，圖4給出了19個(gè)探空站基于兩套模型的預(yù)報(bào)K值與積分法K值(當(dāng)作真值)的時(shí)序?qū)Ρ?

由圖4可知，在低海拔地區(qū)(<600 m，上海、海口、北京、章丘、百色、懷化、漢中)和高海拔地區(qū)(>1500 m，西昌、昆明、西寧、甘孜、拉薩、玉樹、那曲)，兩套模型的預(yù)報(bào)K值相差較大，最大達(dá)到0.004.而在中海拔地區(qū)(600～1500 m，伊寧、太原、烏魯木齊、呼和浩特、貴陽)，兩套模型的預(yù)報(bào)K值偏差很小，尤其是太原、貴陽和呼和浩特等站，均僅相差0.0001，可見在我國(guó)中海拔地區(qū)，高程對(duì)Emardson精化模型影響很小.

圖4 Emardson-I和Emardson-H模型的2019年預(yù)報(bào)K值時(shí)序變化Fig.4 Time series variation of predicted K value based on Emardson-I and Emardson-H models in 2019

就兩套K值預(yù)報(bào)模型的時(shí)序變化而言，積分法K值與Emardson-H模型的一致性更好，并且海拔越高，Emardson-H模型預(yù)報(bào)的K值精度也會(huì)越高.就中海拔地區(qū)而言，Emardson-H模型與Emardson-I模型預(yù)報(bào)K值與積分法K值相差很小.

表4統(tǒng)計(jì)了高中低不同海拔19個(gè)探空站的K值預(yù)報(bào)精度.從表中可知，中國(guó)地區(qū)Emardson-H模型的適用性相對(duì)較好，所得結(jié)論與圖4類似.低海拔地區(qū)的Emardson-H模型的K值預(yù)報(bào)精度較Emardson-I模型略有提高，其預(yù)報(bào)K值的平均絕對(duì)誤差(MAE)均值和均方根(RMS)均值分別為0.00144和0.00173，均優(yōu)于Emardson-I預(yù)報(bào)模型的0.00183和0.00218；在中海拔地區(qū)，Emardson-H模型相對(duì)于Emardson-I模型的K值預(yù)報(bào)精度提升很小，RMS均值從0.00208提高到0.00202；而在高海拔地區(qū)，Emardson-H模型相對(duì)于Emardson-I模型的預(yù)報(bào)精度得到了顯著的提升，RMS均值從0.00312提高到0.00181.從表4中還可以看出，用Emardson-I模型預(yù)報(bào)的K值RMS變化較大，在0.001～0.006之間，海拔越高，偏差越大；而利用Emardson-H模型預(yù)報(bào)的K值RMS集中在0.002左右，穩(wěn)定性比前者有大幅度改善.

2.3 GNSS-PWV反演精度分析

2.3.1 GNSS-PWV時(shí)序變化及其偏差

由表1和表2可知，安慶、銀川和拉薩等處的探空站和GNSS站距離較近，高程相差最多38 m，可認(rèn)為探空站與GNSS站所受到的對(duì)流層延遲非常接近，可用探空PWV作為真值來檢驗(yàn)基于Emardson-I和Emardson-H模型的GNSS-PWV精度和可靠性.

限于篇幅，本文僅選擇6個(gè)不同海拔的探空站和GNSS并址站，分析利用Emardson-I和Emardson-H模型的預(yù)報(bào)K值反演的2019年GNSS-PWV和探空PWV時(shí)序變化.如圖5所示，兩組GNSS-PWV與探空PWV在變化趨勢(shì)上表現(xiàn)一致，PWV波峰出現(xiàn)在夏季，PWV波谷出現(xiàn)在冬季，峰值和谷值基本吻合，符合中國(guó)夏季多雨、冬季干燥的氣候特點(diǎn)，說明兩套模型反演的GNSS-PWV具有較高的一致性.

由兩套模型計(jì)算的GNSS-PWV與探空PWV之間的偏差圖可看出，基于Emardson-H模型的GNSS-PWV偏差均小于Emardson-I模型，且夏季的GNSS-PWV偏差較大.在高海拔地區(qū)(昌都、拉薩)，GNSS-PWV偏差范圍相對(duì)較小(-0.5～0.5 mm)；在中海拔地區(qū)(銀川、若羌)，基于Emardson-I和Emardson-H模型的GNSS-PWV偏差范圍相對(duì)較大(-1.0～0.5 mm)；而在低海拔地區(qū)(安慶、青島)，PWV偏差范圍相對(duì)更大(-2.0～0.5 mm)，但基于Emardson-H模型的GNSS-PWV更接近探空PWV，比Emardson-I模型有更為明顯的改善.

表4 2019年基于Emardson-I和Emardson-H模型的K值預(yù)報(bào)精度Table 4 Prediction accuracy of K value based on Emardson-I and Emardson-H models in 2019

2.3.2 GNSS-PWV精度

由表5可以得到，低海拔地區(qū)(青島、安慶)基于Emardson-I模型的GNSS-PWV的MAE均值為0.629 mm，RMS均值為0.752 mm，而基于Emardson-H模型的GNSS-PWV的MAE和RMS均值分別為0.433 mm和0.530 mm，RMS降低了0.222 mm；中海拔地區(qū)(若羌、銀川、貴陽)基于Emardson-I模型的GNSS-PWV的MAE均值為0.274 mm，RMS均值為0.342 mm，而基于Emardson-H模型的GNSS-PWV的MAE和RMS均值分別為0.266 mm和0.333 mm，RMS僅降低了0.009 mm；而高海拔地區(qū)(和田、都蘭、昌都、甘孜、拉薩)基于Emardson-I模型的GNSS-PWV的MAE均值為0.159 mm，RMS均值為0.213 mm，而基于Emardson-H模型的GNSS-PWV的MAE和RMS均值分別為0.119 mm和0.154 mm，兩種模型的RMS相差在0.1 mm以內(nèi).

綜上所述，在中國(guó)地區(qū)，基于Emardson-H模型的GNSS-PWV反演精度明顯高于Emardson-I模型，其MAE和RMS均值均優(yōu)于Emardson-I模型.在中高海拔地區(qū)，基于Emardson-H模型的GNSS-PWV相對(duì)于Emardson-I模型，精度有所提高，但提高有限；相比而言，Emardson-H模型在低海拔地區(qū)的改正效果更為明顯，這是因?yàn)榈秃０蔚貐^(qū)的ZTD和PWV值均較大，高精度K值更加有利于提高GNSS-PWV精度，說明在低海拔地區(qū)考慮高程因子同樣重要.

3 討論與結(jié)論

水汽作為重要的氣象要素之一，是災(zāi)害性天氣形成的重要因子，在降水監(jiān)測(cè)和天氣預(yù)報(bào)中扮演著關(guān)鍵角色.中國(guó)幅員遼闊、地形起伏大、探空站時(shí)空分辨率低，導(dǎo)致大部分區(qū)域的水汽資料比較匱乏.地基GNSS技術(shù)能探測(cè)高時(shí)空分辨率PWV，而K值是從GNSS-ZTD轉(zhuǎn)換到GNSS-PWV的重要參數(shù).因此，本文基于Emardson模型建立了兩套適合中國(guó)地區(qū)的高精度K值模型，重要結(jié)論如下：

(續(xù)圖5)

圖5 2019年基于Emardson-I和Emardson-H模型的GNSS-PWV時(shí)序變化及其偏差Fig.5 Time series variations and their deviations of GNSS-PWV based on Emardson-I and Emardson-H models in 2019

表5 2019年基于Emardson-I和Emardson-H模型的GNSS-PWV精度統(tǒng)計(jì)Table 5 The accuracy statistics of GNSS-PWV based on Emardson-I and Emardson-H models in 2019

(1)K值建模需要考慮緯度和高程這兩個(gè)關(guān)鍵參數(shù).K值與測(cè)站緯度和高程之間呈負(fù)相關(guān)，其相關(guān)系數(shù)分別為0.735和0.941，而與測(cè)站經(jīng)度的相關(guān)性較小.

(2)Emardson-H模型的K值預(yù)報(bào)精度優(yōu)于Emardson-I模型，在高海拔地區(qū)尤為明顯.Emardson-H模型預(yù)報(bào)的K值MAE和RMS均值分別為0.00150和0.00182，均優(yōu)于Emardson-I模型的0.00221和0.00255.

(3)基于Emardson-H模型的GNSS-PWV精度也優(yōu)于Emardson-I模型，尤其在低海拔和高海拔地區(qū).基于Emardson-H模型的GNSS-PWV的MAE和RMS均值分別為0.226 mm和0.283 mm，均優(yōu)于基于Emardson-I模型的0.288 mm和0.360 mm.

綜上所述，就K值而言，Emardson-H模型預(yù)報(bào)的K值精度優(yōu)于Emardson-I模型，尤其在高海拔地區(qū)；而就GNSS-PWV來說，不僅在高海拔地區(qū)，并且在低海拔地區(qū)，Emardson-H模型反演的GNSS-PWV精度均優(yōu)于Emardson-I模型.這是因?yàn)榈秃０蔚貐^(qū)的ZTD值較大，高精度K值更有利于提高PWV反演精度，說明本文所建立的全國(guó)K值Emardson模型適用于大部分地區(qū).但分析圖4可以發(fā)現(xiàn)，某些探空站(比如海口、伊寧、烏魯木齊、西昌等)在春冬季節(jié)的模型預(yù)報(bào)K值偏差較大，意味著K值模型不僅僅要考慮年積日、緯度和高程，可能也要考慮當(dāng)?shù)氐募竟?jié)變化和氣候特征；此外，本文僅考慮了影響因子的線性特征，而沒有考慮各參數(shù)的非線性變化特點(diǎn)，后續(xù)我們還可基于以上兩點(diǎn)進(jìn)一步完善該模型.

致謝感謝美國(guó)懷俄明大學(xué)天氣數(shù)據(jù)網(wǎng)和中國(guó)地震局GNSS產(chǎn)品服務(wù)平臺(tái)為本文提供了探空站資料和GNSS數(shù)據(jù).