李得宴,楊維芳,高志鈺,4,李蓉蓉

(1.蘭州交通大學 測繪與地理信息學院,蘭州 730070;2. 地理國情監(jiān)測技術應用國家地方聯(lián)合工程研究中心,蘭州730070;3. 甘肅省地理國情監(jiān)測工程實驗室,蘭州730070;4. 中國地震局地質研究所,北京 100029)

0 引 言

水汽作為氣象學和天氣預報中關注的一個重要因素,其變化驅動著天氣、氣候產(chǎn)生相應的變化,同時也參與了地球氣候系統(tǒng)的能量流動和水循環(huán)[1]．Bevis等[2]在1992年提出了基于GPS的大氣水汽探測方法,為大氣水汽探測提供了一種新方法．大氣加權平均溫度Tm作為地基全球衛(wèi)星導航系統(tǒng)(GNSS)水汽反演過程中的一個重要變量,可以利用探空站資料通過數(shù)值積分的方法精確獲得．但是計算Tm需要飽和水汽壓Es參與其中,氣象學研究中各研究人員建立了多種飽和水汽壓模型,目前在國內(nèi)地基GNSS水汽反演研究中計算飽和水汽壓主要選取了其中的三種模型(Magnus-Tetens模型、BUCK模型、Goff-Gratch模型),查閱資料發(fā)現(xiàn)這三種飽和水汽壓模型無論是結構還是計算結果都存在一定的差異,可能會對大氣可降水量(PWV)的準確反演產(chǎn)生影響, 這將不利于不同研究人員研究結果的對比和氣象數(shù)據(jù)的融合．所以本文以香港為研究區(qū),利用上述三種不同的飽和水汽壓模型計算飽和水汽壓Es,通過數(shù)值積分的方法計算Tm;利用GAMIT解算得到了旱雨兩季(2、7月)的天頂濕延遲(ZWD),然后利用Matlab編程模擬了由ZWD轉化為PWV的過程．代入探空站數(shù)值積分求得的Tm和GAMIT解算得到的ZWD并最終計算得到每天的PWV,進一步對比三種不同飽和水汽壓模型參與計算得到的Tm和PWV,最后分析這三種飽和水汽壓模型是否均適用于地基GNSS水汽反演．由于在研究過程中發(fā)現(xiàn)部分研究人員將Buck模型中的變量T作為露點溫度來進行相關研究,所以本文就這一問題進行了探究,并闡述了該做法的不合理性,這將為地基GNSS水汽反演中的部分研究提供參考．

1 GNSS水汽反演原理

1.1 PWV的計算步驟

在GNSS氣象學中,利用GAMIT軟件解算GNSS觀測數(shù)據(jù)可以得到對流層天頂總延遲(ZTD),利用干延遲模型可以計算得到天頂對流層干延遲(ZHD),ZTD減去ZHD就可以得到對流層天頂濕延遲(ZWD)．

ZWD=ZTD-ZHD.

(1)

最后通過轉換系數(shù)Π將ZWD轉換為PWV．

PWV=Π×ZWD.

(2)

轉換系數(shù)Π為

(3)

1.2 大氣加權平均溫度的計算

大氣加權平均溫度Tm的計算方法有以下三種方法:①常數(shù)法；②探空資料數(shù)值積分法；③模型法.本文將選擇不同飽和水汽壓模型在氣溫t下計算得到有差異的飽和水汽壓值并通過數(shù)值積分(方法②)來計算得到不同的Tm，然后進行進一步計算、對比和分析。探空資料數(shù)值積分計算Tm公式為

(4)

式中，ei為傳感器隨探空氣球上升過程中測得的第i和i+1層水汽壓的平均值,但是水汽壓e的值不能直接由傳感器測量得到,需要通過測得的相對濕度RH和測得的第i和i+1層的平均溫度Ti所對應的飽和水汽壓值Es計算得到(如式5)．

(5)

式中：e為飽和水汽壓；RH為相對濕度;Es為飽和水汽壓.

2 飽和水汽壓計算模型

飽和水汽壓是定量的空氣在給定的溫度下所能保持的最大的水汽壓,溫度越高,空氣達到飽和所需的水汽就越多．當水汽壓高于飽和水汽壓時,水汽將會凝結為小水滴．下面的三種模型是目前國內(nèi)研究人員在地基GNSS水汽反演的研究中所采用的不同的飽和水汽壓模型．

2.1 Magnus-Tetens模型

(6)

(7)

式(6)為相對于純水面的飽和水汽壓公式(適用溫度為-45 ℃～60 ℃),一般情況下,水會在 0 ℃以下結晶成冰,但是研究發(fā)現(xiàn)純水可以在-48 ℃時還以液態(tài)存在而形成過冷水,并且在層云和積云中存在過冷水,所以式(6)雖然為相對于水面的飽和水汽壓公式,但將其溫度適用范圍擴展到0 ℃以下,使其也可用于過冷水狀態(tài)下的飽和水汽壓的計算[3-5]．式(7)為相對于純冰面的飽和水汽壓公式(適用溫度為-60 ℃～0 ℃),式(6)、式(7)中,T為大氣溫度(℃),Esm代表該模型在相應溫度下求得的飽和水汽壓值(hPa),該模型由德國氣象學家O. Tetens在1930年提出[6]．1967年Murray對其進行精度的評定,發(fā)現(xiàn)該公式精度較高,可以被應用于氣象學中,但是其相對精度會隨著溫度的降低而降低[7]．文獻[8-9]在地基GNSS研究的過程中均應用了該公式,但是在應用過程中并沒有根據(jù)溫度明確地區(qū)分應用水面、冰面公式,而是只用了相對于水面的飽和水汽壓公式來計算探空氣球上升過程中的每一層的飽和水汽壓．

2.2 Goff-Gratch模型

(8)

(9)

式(8)為相對于純水面飽和水汽壓公式(適用溫度為0 ℃～100 ℃),式(9)為相對于純冰面的飽和水汽壓公式(適用溫度為-100 ℃～0 ℃),式(8)、(9)中,Esg代表該模型求得的飽和水汽壓值(hPa),T為大氣溫度(K),T1為水的三相點溫度(273.16 K)．該模型由Goff和Gratch于1946年提出,Goff在1957年對該公式做了修訂后被世界氣象組織(WMO)作為計算飽和水汽壓的推薦公式[10]．式(8)的適用溫度為(0 ℃～100 ℃),但在WMO(2000)的技術規(guī)范中注明該模型在溫度為-50 ℃的過冷水面上仍然可以使用,且誤差極小[4-5,11]．但是由于該模型建立之初并未考慮過冷水狀態(tài),所以本文應用模型時該時仍以0 ℃作為分界點．文獻[11-12]在地基GNSS水汽反演的研究中利用該模型來計算飽和水汽壓,但僅應用了相對于水面的飽和水汽壓模型來參與Tm的計算[12-13]．

2.3 BUCK模型

(10)

(11)

式(10)為相對于純水面的飽和水汽壓公式(適用溫度為-45 ℃～60 ℃),所以式(10)也適用于過冷水狀態(tài)下的飽和水汽壓的計算,式(11)為相對于純冰面的飽和水汽壓公式(適用溫度為-65 ℃～0 ℃),式(9)、(10)中,Esb代表該模型求得的飽和水汽壓值(hPa);T為大氣溫度(℃)．該模型由Buck于1981年利用極小極大曲線擬合(Minimax fitting)法將不同溫度區(qū)間的最大擬合誤差最小化,進而得到該模型的最優(yōu)系數(shù)[4-5,14]．文獻[15-18]在利用數(shù)值積分法計算Tm時均使用該模型來得到所需的飽和水汽壓值[15],但是目前在地基GNSS反演中應用該模型時,有大部分研究人員將T當作露點溫度Td(℃)來計算[16-18],所以下面就該模型分別將T作為大氣溫度和露點溫度的兩種情況進行討論,并對比這樣的處理方式對Tm產(chǎn)生的影響,下文中把T當作露點溫度來處理的BUCK模型稱為BUCK 2模型．

3 研究區(qū)域及其方法

本文研究區(qū)域為香港,地處亞熱帶季風氣候區(qū),分旱雨兩季,5—10月降雨較多,其余月份降雨較少[19]．香港king’s park探空站與香港地區(qū)衛(wèi)星定位參考網(wǎng)中的昂船洲(HKSC)站相距約3 km,距離較近;海拔高度相差約40 m,地形起伏差異較?。话闱闆r下,利用GAMIT解算PWV時所用到的Tm都是通過建立的Tm模型來得到,但是為了避免建立模型環(huán)節(jié)帶來的誤差導致直接觀察到不同飽和水汽壓公式參與計算引起的PWV差異,本文近似地使用king’s park探空站的2016年旱(2月)雨(7月)兩季高空氣象數(shù)據(jù)計算得到的Tm來代替HKSC站進行解算時所用到的Tm．選取第二小節(jié)所述的不同飽和水汽壓模型計算得到的飽和水汽壓Es參與Tm的計算,將計算得到的不同的Tm分別應用到HKSC站點的GNSS水汽反演中并得到不同的大氣可降水量PWV值,將king’s park探空站的PWV視為真值,對比在反演過程中不同飽和水汽壓模型值參與計算得到的PWV值的差異,并分析差異存在的原因．探空站計算大氣可降水量的公式如式(12)所示,式(13)為式(12)的數(shù)值積分式:

(12)

(13)

(14)

式中:g(m/s2)為地球重力加速度，m/s2;pi為不同高度觀測層的氣壓,hPa;qi為不同高度觀測層的比濕,g/kg,其計算公式為式(14);e為水汽壓．

分別查閱香港天文臺2013、2016版的香港氣象及潮水觀測摘要了解到,隸屬于香港天文臺的king’s park無線電探空站所使用的無線電探空儀是由芬蘭Vaisala公司研發(fā)制造的．且該探空站在2006年7月之前使用Vaisala WS80型探空儀,2006年7月1日至2016年11月使用Vaisala RS92型探空儀,2016年11月至今使用Vaisala RS41型探空儀[20]．查閱Vaisala公司關于探空儀的技術支持報告,了解到Vaisala系統(tǒng)在計算飽和水汽壓時采用了國際水和水蒸氣性質協(xié)會 (IAPWS)于1995年推薦的飽和水汽壓公式,如式(15)、(17)所示[21],式(15)為相對于水面的飽和水汽壓公式(適用溫度為0 ℃～373 ℃),式(17)為相對于冰面的飽和水汽壓公式(適用溫度為-100 ℃～0.01 ℃)．式(15)由A.Saul和W.Wagner于1987年提出,W.Wagner于1992年參照ITS-90(1990國際溫標)進行了修訂得到了式(15)[22],以下將式(15)、(17)簡稱為W.Wagner模型．

c4v3.5+c5v4+c6v7.5),

(15)

(16)

(17)

(18)

式(15)中：es為飽和水汽壓;pc=22.064 hPa,為零界壓力;Tc=647.097 K,為零界溫度;T為大氣溫度, K; 系數(shù)c1=-7.859 517 83;c2=1.844 082 59;c3=-11.786 649 7;c4=22.680 741 1;c5=-15.961 871 9;c6=1.801 225 02;v的計算如式(16)所示．式(17)中,pn=6.116 57 hPa,為三相點溫度的水汽壓;系數(shù)a0=-13.928 169;a1=34.707 823,θ的計算如式(18)所示,其中T為大氣溫度,K;Tn=273.16 K,為水的三相點溫度．

雖然在地基GNSS水汽反演的各研究中并未將式(15)、(17)應用到飽和水汽壓的計算過程中,但是由于在對比時將香港king’s park探空站的PWV作為真值來參考,由上所述可知式(15)、(17)參與了探空站的PWV的計算,所以在下面的討論中式(15)、(17)也將參與地基GNSS水汽反演過程中的Tm的計算,并和其他模型進行比較．

此次反演對比試驗將會分別選擇旱雨兩季降水量差異較大的2月、7月份分別進行,且在基線解算時引入LHAZ、PIMO、SHAO、CUSV四個IGS站．IGS站、探空站、衛(wèi)星定位參考網(wǎng)中的HKSC站的地理位置分布如圖1所示．

圖1 IGS站、CORS站、king’s park 探空站分布

4 不同飽和水汽壓模型對Tm及PWV的影響

因為選擇不同的飽和水汽壓模型MT模型、GG模型、BUCK模型、BUCK 2模型、WW模型會首先影響大氣加權平均溫度Tm(下面的分析中將各不同飽和水汽壓模型參與計算得到的Tm簡寫為Tm(MT)、Tm(GG)、Tm(BUCK)、Tm(BUCK 2)、Tm(WW)),進而影響PWV的反演精度,所以在討論不同飽和水汽壓模型對PWV精度的影響之前會先對比利用不同飽和水汽壓模型計算得到的Tm的差異,再將GNSS最終反演得到的PWV值與探空站PWV值對比分析．

4.1 飽和水汽壓模型的選取對大氣加權平均溫度(Tm)的影響

由于大氣水汽幾乎全部集中在氣象學定義的對流層(0～12 km)內(nèi),所以在利用數(shù)值積分法計算Tm時僅采用探空數(shù)據(jù)中的大氣高度為0～12 km記錄的數(shù)據(jù)．通過分析香港地區(qū)2016年低溫月份(2月)、高溫月份(7月)在最接近大氣層12 km處的溫度發(fā)現(xiàn):2月大氣層頂最低溫度為-56.5 ℃,最高溫度為-44.5 ℃,平均溫度為-49.9 ℃;7月大氣層頂最低溫度為-52.9 ℃,最高溫度為-43.5 ℃,平均溫度為-47.19 ℃．目前的地基GNSS水汽反演的研究中所采用的飽和水汽壓模型以及香港探空站使用的Vaisala探空儀采用的W.Wagner模型已分別在第2、3小節(jié)中做了簡略介紹,且由上述介紹可知Magnus-Tetens模型、BUCK模型的水面和冰面公式的適用溫度范圍存在重疊,文獻[12]就存在以0 ℃作為冰、水面公式的分界點,還是以-45 ℃作為冰、水面公式的分界點的問題進行了實驗驗證,實驗表明:利用探空儀測定的相對濕度計算水汽壓時誤用冰面飽和水汽壓公式對對流層層頂(12 km左右)的影響在高精度需求中是不容忽視的[23]．世界氣象組織(WMO)建議在0 ℃以下的相對濕度的測定應相對于水面,且Vaiaala探空儀在0 ℃以下的相對濕度的測定也是相對于水面的[24]．綜上原因,本文在計算飽和水汽壓時選擇-45 ℃作為Magnus-Tetens模型、BUCK模型的水面和冰面公式的溫度分界點．圖2是由不同飽和水汽壓模型計算得到2月、7月的Tm變化趨勢圖．對比圖2中的(a)、(b)兩張圖可以看出,7月的Tm明顯高于2月．單獨分析圖(a)、(b),可以發(fā)現(xiàn)在其余氣象因子都相同的情況下,利用不同飽和水汽壓模型通過式(4)數(shù)值積分計算得到的Tm變化趨勢相同;但是由BUCK 2模型(以露點溫度Td作為變量的飽和水汽壓模型)參與計算得到的Tm與其他飽和水汽壓模型參與計算得到的Tm相差較大,并且七月份的差異大于二月份．

(a)2月份 (b)7月份圖2 利用不同飽和水汽壓模型通過數(shù)值計算法得到的香港2、7月份Tm變化趨勢

各飽和水汽壓參與計算得到的Tm月平均值m如式(19)所示,計算結果如表1所示．

(19)

式中：m為單個飽和水汽壓模型參與Tm計算得到的Tm的平均值；Tmi為探空氣球第i次升空探測并計算得到的Tm值；n為探空氣球升空探測次數(shù)．

表1 不同飽和水汽壓模型參與計算得到的香港2、7月份Tm的均值m K

通過圖2中的(a)、(b)兩張圖及相應的表1可以看出Tm(BUCK)、Tm(MT)、Tm(GG)、Tm(WW)相差不大,一致性較好,但是Tm(BUCK 2)與其他的Tm值相差較大,存在明顯的系統(tǒng)性差異．下面將結合圖3所示的飽和曲線圖分析這種較為明顯的差異性的存在原因．圖3的橫軸為大氣溫度T,縱軸為飽和水汽壓Es,紅色曲線為飽和曲線,在飽和曲線上的水汽處于飽和狀態(tài),水的蒸發(fā)速率等于凝結速率,水在氣態(tài)、液態(tài)兩種狀態(tài)之間的轉換處于動態(tài)平衡．在曲線上的相對濕度為100%,所以此時的飽和水汽壓等于實際水汽壓(Es=e),且大氣溫度等于露點溫度(T=Td)．相對濕度、飽和水汽壓、水汽壓的關系為

(20)

在曲線下方的區(qū)域A中,水汽處于不飽和狀態(tài),水的蒸發(fā)速率大于凝結速率,相對濕度小于100%,如果水汽壓e不變,而想要處于A區(qū)域的不飽和狀態(tài)3變?yōu)轱柡蜖顟B(tài)2,就需要將狀態(tài)3的溫度(氣溫T)降低到狀態(tài)2的溫度(露點溫度Td),使得飽和水汽壓值Es的值降低,式(22)中的相對濕度RH就會達到100%,既水汽達到飽和狀態(tài)．在曲線上方的區(qū)域B中,水汽處于過飽和狀態(tài),水的凝結速率大于蒸發(fā)速率,如果水汽壓e不變,使處于B區(qū)域的過飽和狀態(tài)1變?yōu)轱柡蜖顟B(tài)2,就需要將狀態(tài)1的溫度(氣溫T)升高到狀態(tài)2的溫度(露點溫度Td)．然而分析發(fā)現(xiàn)由king’s park探空站的探空儀所測得的相對濕度數(shù)據(jù)中只有個別值能達到100%,且相對濕度隨探空儀高度的上升有下降趨勢,所以由上述分析可知,探空儀上升過程所觀測的水汽并沒有全部達到飽和狀態(tài),而用露點溫度Td作為BUCK飽和水汽壓公式的變量,得到的就是水汽為處于飽和狀態(tài)時的水汽壓,就會高于實際值．所以就會導致Tm產(chǎn)生如上所述的明顯差異,使得無論是Tm(BUCK 2)的平均值還是標準差都高于其他飽和水汽壓模型參與計算得到的Tm．

圖3 飽和曲線

4.2 飽和水汽壓模型的選取對PWV的影響

因為在GAMIT中為了一般化和方便大量計算,通常是利用文件中默認的以測站溫度為變量的Bevis線性Tm模型或者是將其替換為較精確的本地化Tm模型來計算Tm,但為了避免建立模型環(huán)節(jié)帶來的誤差而導致直接觀察不同飽和水汽壓公式計算得到的值參與反演引起的PWV的差異,本文直接將距離HKSC較近的king’s park探空站所求得的Tm值作為計算轉換系數(shù)Π時的輸入值,只利用了GAMIT解算得到的濕延遲ZWD,將反演后續(xù)的轉換系數(shù)Π及PWV的計算過程利用Matlab編程實現(xiàn)．且由上述分析,采用BUCK 2模型是不合理的,會引起較大的誤差,所以在下面的PWV的對比分析中將不再討論BUCK 2模型參與計算引起的PWV的差異．不同飽和水汽壓模型參與計算得到的PWV的趨勢圖如圖4所示(下面的分析中將各不同飽和水汽壓模型參與計算得到的PWV簡寫(PWVMT、PWVGG、PWVBUCK、PWVWW):

由圖4中的(a)、(b)圖可以看出,GNSS反演得到的PWV的趨勢與Radiosond探空站數(shù)據(jù)中PWV趨勢一致,不同飽和水汽壓參與計算得到的PWV差異較?。蕴娇照綪WV數(shù)據(jù)為真值,計算了PWVMT、PWVGG、PWVBUCK、PWVWW的平均偏差(MEAN)、平均絕對誤差(MAE)、均方根誤差(RMSE),結果如表2、3所示.

(a)2月份 (b)7月份圖4 2、7月份大氣可降水量(PWV)趨勢圖

表2 2月份PWV精度驗證 mm

表3 7月份PWV精度驗證 mm

進一步通過獨立樣本檢驗來研究討論四種不同的飽和水汽壓模型參與計算得到的PWV與探空站真值的差異性．在進行數(shù)據(jù)差異性檢驗之前,應先通過W檢驗(Shapio-Wik檢驗)來檢驗各組數(shù)據(jù)是否近似服從正態(tài)分布,然后再決定使用哪一種獨立樣本檢驗．當一組數(shù)據(jù)的W檢驗的顯著性P值小于0.05時,否定原假設,認為該組數(shù)據(jù)不近似服從正態(tài)分布．結果表明,除7月的探空站PWV的W檢驗的P值為0.2明顯大于0.05以外,其余各組數(shù)據(jù)的W檢驗的P值均小于0.05．說明除了7月的PWV(Radiosond)近似服從正態(tài)分布,2月、7月的其他組的PWV均不知道其分布狀態(tài),又因為Mann Whitney檢驗屬于非參數(shù)檢驗,進行該檢驗前并不需要知道總體的分部情況,所以本文將采用Mann Whitney檢驗來判斷兩組獨立的數(shù)據(jù)之間是否存在顯著差異．該檢驗通過檢驗兩組數(shù)據(jù)的均值是否存在差異來判斷兩個獨立的樣本是否存在顯著差異,即若兩個樣本有差異,則他們的中心位置就會不同．當一組GNSS反演得到的PWV值與探空站的PWV值的Mann Whitney檢驗的顯著性P值大于0.05時,接受原假設,認為反演得到的PWV與真值(探空站PWV值)不存在明顯差異．2月、7月的GNSS反演PWV值與真值的Mann Whitney檢驗的顯著性P值分別列在表4中.

表4 Mann Whitney獨立樣本檢驗的顯著性P值

由表3、4可以看出,2、7月的PWV反演精度都較高,且整體上2月份的PWV反演精度要高于7月份．但是反演得到的PWV的精度差異不大．在表5中,進一步通過Mann Whitney檢驗的顯著性P值可以看出,無論是2月還是7月,顯著性P值均大于0.05,所以可以進一步明確各不同飽和水汽壓參與計算得到的PWV與真值之間不存在具有統(tǒng)計意義的顯著性差異．這說明各飽和水汽壓公式雖然存在形式和計算結果上的差異,但是該差異并沒有使得GNSS反演得到的PWV值產(chǎn)生明顯的不同．

5 結 論

本文針對目前地基GNSS水汽反演研究當中計算Tm的過程中所應用的三種不同的飽和水汽壓模型,對比在反演過程中不同飽和水汽壓模型參與計算得到的飽和水汽壓值所引起的中間變量Tm和最終反演得到的PWV值的差異,評價這三種飽和水汽壓模型能否較好地應用于地基GNSS水汽反演的研究中,得到如下的結論:1)Magnus-Tetens模型、Buck模型、Goff-Gratch模型應用在地基GNSS水汽反演過程中都能得到較好的反演精度,且利用三個不同模型反演得到的PWV值的差異性也較?。肰aiaala探空系統(tǒng)采用的飽和水汽壓模型參與GNSS水汽反演得到的PWV值和中間值Tm與上述三種模型參與計算得到的值有很好的一致性．所以利用數(shù)值積分求Tm值并用來建立Tm模型時,這三個不同的飽和水汽壓模型都可以被用來提供計算樣本Tm時所用到的飽和水汽壓．2)在使用BUCK模型時,誤將T作為露點溫度來計算飽和水汽壓是不科學、不合理的,將會產(chǎn)生較大的誤差,所以以后的研究中在使用該公式計算飽和水汽壓時,應該將T當作大氣溫度,才能利用BUCK模型得到合理的飽和水汽壓值,進而應用到GNSS水汽反演過程中．