莫智翔 李佳豪 周呂 黃良珂 劉立龍

0 引言

水汽及其變化是天氣、氣候變化的主要驅(qū)動力,是災(zāi)害性天氣形成和演變過程中的重要因子．大氣水汽的變化與降水直接相關(guān),在大氣能量傳輸、天氣系統(tǒng)演變、大氣輻射收支、全球氣候變化等多種氣象演變中扮演著重要的角色[1-2]．大氣可降水量(Precipitable Water Vapor,PWV)又稱大氣水汽總量,定義為地面上大氣柱中的總水汽量,是表征大氣水汽含量以及空中水資源的重要指標(biāo)[3-4]．目前,對大氣水汽含量的探測主要依靠常規(guī)的無線電探空站、雷達觀測和水汽輻射計等技術(shù),但使用費用昂貴,且時空分辨率低,離監(jiān)測和預(yù)報中小尺度災(zāi)害性天氣的要求還有很大差距．全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System,GNSS)能較好地彌補傳統(tǒng)大氣水汽探測技術(shù)在時空分辨率上的不足,并提供精細化氣象預(yù)報所需要的高精度、大容量、近實時的大氣水汽資料．

GNSS反演的PWV的計算通常由GNSS衛(wèi)星信號傳播延遲時的對流層天頂濕延遲(Zenith Wet Delay,ZWD)乘以水汽轉(zhuǎn)換系數(shù)得到,而水汽轉(zhuǎn)換系數(shù)的關(guān)鍵參數(shù)就是大氣加權(quán)平均溫度(Tm)．因此,在GNSS反演大氣水汽中,Tm是獲取GNSS大氣水汽的一個關(guān)鍵參數(shù),準(zhǔn)確地計算Tm值是提高PWV計算精度、實現(xiàn)實時PWV估計的重要保證．Tm值可以利用探空資料或大氣再分析產(chǎn)品等氣象參數(shù)通過數(shù)值積分法計算獲取,但受限于地理、數(shù)據(jù)獲取等因素,不容易獲取任意位置處的Tm值．為此通常需要建立一個準(zhǔn)確的Tm模型來滿足這些需求,為用戶提供便利．目前,Tm模型主要分為需要實測氣象參數(shù)的Tm模型和不需要實測氣象參數(shù)的Tm模型．在需要氣象參數(shù)的Tm模型中,文獻[5]利用北美8 718個探空站發(fā)現(xiàn)Tm與地表溫度(Ts)具有良好的線性關(guān)系,并提出了一種常用的Tm計算經(jīng)驗?zāi)Ｐ?Bevis模型)．由于Bevis模型應(yīng)用到其他地區(qū)會存在一定的系統(tǒng)偏差[6],諸多學(xué)者針對不同地區(qū)對Bevis模型進行分析和改進[7-8]．文獻[9]利用探空數(shù)據(jù)建立了適用于中國東部地區(qū)的Tm模型,并取得良好的效果;文獻[10]分析了Tm與Ts、緯度、海拔、水汽壓和大氣壓之間的關(guān)系,建立了Tm單因子和多因子回歸模型；文獻[11]則在Bevis公式的基礎(chǔ)上在中國西部地區(qū)建立了一種顧及Ts、高程和季節(jié)變化的Tm模型．在不需要氣象參數(shù)的Tm模型中,這類模型一般是基于局部或全球多年的Tm數(shù)據(jù)構(gòu)建的經(jīng)驗?zāi)Ｐ停墨I[12]首先構(gòu)建了歐洲地區(qū)與緯度和季節(jié)相關(guān)的非氣象參數(shù)Tm模型(Emardson模型),諸多學(xué)者也對其進行了改進[13]．文獻[14]提出一個全球經(jīng)驗對流層格網(wǎng)模型(GPT2w模型),該模型有1°×1°(GPT2w-1)和5°×5°(GPT2w-5)兩種分辨率,只要輸入目標(biāo)經(jīng)緯度、高程和相應(yīng)的年積日,即可提供所需的包括Tm在內(nèi)的對流層參數(shù)．文獻[15]建立了中國區(qū)域顧及垂直遞減率函數(shù)的非氣象參數(shù)Tm模型．雖然非氣象參數(shù)Tm模型不用依賴實測的氣象參數(shù),但與采用實測氣象參數(shù)的Tm模型相比,存在一定的系統(tǒng)誤差,其精度在局部地區(qū)的表現(xiàn)還有待進一步驗證．

青藏高原位于歐亞大陸中部,是世界海拔最高的高原,被稱為“世界屋脊”．它的冷熱作用對中國、東亞乃至世界的大氣環(huán)流和氣候變化具有非常重要的影響[16]．青藏高原地形多樣、水汽分布復(fù)雜,且探空站分布稀疏,導(dǎo)致青藏高原水汽的監(jiān)測成為大氣探測中的難題．此外,大多現(xiàn)有的Tm模型在青藏高原地區(qū)的適用性較差,且未在該地區(qū)有深入的研究．因此,研究建立高精度的青藏高原區(qū)域Tm經(jīng)驗?zāi)Ｐ途哂兄匾默F(xiàn)實意義．本文利用2014—2017年分布在青藏高原地區(qū)的探空資料分析了Tm與Ts、緯度和高程在青藏高原地區(qū)的相關(guān)性,在此基礎(chǔ)上建立了適用于青藏高原地區(qū)的顧及時空變化因素的多因子Tm模型(TPTm模型),并利用2018年的探空數(shù)據(jù)作為參考,驗證該模型的精度,為青藏高原地區(qū)高精度的GNSS水汽反演、應(yīng)用提供重要的參考依據(jù)．

1 計算原理及數(shù)據(jù)來源

1.1 計算原理

GNSS反演的大氣水汽PWV[5]由天頂濕延遲ZWD乘以水汽轉(zhuǎn)換系數(shù)Π得到:

VPW=Π·ZWD，

(1)

(2)

式中:ρw=1×103kg/m3為液態(tài)水的密度;Rv=461.495 J·kg-1·K-1為水汽氣體常數(shù);k′2，k3為大氣物理參數(shù),經(jīng)驗值通常分別為22.13±2.20、(3.739±0.012)×105K/hPa;Tm為大氣加權(quán)平均溫度,可以由測站上空水汽壓和絕對溫度沿天頂方向的積分函數(shù)求得,其數(shù)學(xué)表達式[17]為

(3)

式中:e和T分別為測站天頂方向某高度的水汽壓(hPa)和絕對溫度(K);z表示測站上空的垂直高度(m);hs為測站大地高(m)．在探空資料中,一般提供了相對濕度HR和絕對溫度T,而水汽壓e并沒有直接提供,但可以通過飽和水汽壓es(hPa)和露點溫度Td(單位為℃,T=Td+273.15)來計算得到,公式為

(4)

(5)

在實際中,式(3)采用以下積分公式離散化得到:

(6)

Tm=a+bTs,

(7)

式中,公式系數(shù)a和b可利用最小二乘法計算得到．Bevis等[5]利用北美地區(qū)探空資料構(gòu)建的Bevis模型曾被國內(nèi)外諸多學(xué)者采用,其求取的a,b系數(shù)分別取值為70.2和0.72．

1.2 數(shù)據(jù)源介紹

美國懷俄明州立大學(xué)網(wǎng)(http:∥www.weather.uwyo.edu/upperair/sounding.html)可以直接下載探空站的觀測資料．其中探空站資料包括探空氣球每天在0時和12時探測兩次的大氣分層數(shù)據(jù),數(shù)據(jù)包含這兩個時刻的氣壓、溫度、露點溫度、相對濕度等相關(guān)氣象數(shù)據(jù),還提供了測站的位置信息．根據(jù)不同的等壓面,可由數(shù)值積分法計算Tm值,并從探空數(shù)據(jù)中獲得Ts值．本文選取2014—2017年均勻分布在青藏高原地區(qū)13個探空站的觀測數(shù)據(jù)來進行模型構(gòu)建,用2018年的相應(yīng)探空資料進行評估、驗證．探空站點位分布及各站點相關(guān)信息如圖1和表1所示．

圖1 青藏高原地區(qū)13個探空站分布Fig.1 Distribution of the 13 radiosonde stations in Tibetan Plateau

表1 青藏高原地區(qū)13個探空站信息情況

圖2 青藏高原地區(qū)Tm與Ts、高度和緯度的相關(guān)性Fig.2 Relationship between Tm,Ts,altitude and latitude for Tibetan Plateau

2 模型的建立與評估

2.1 加權(quán)平均溫度與地面溫度、高度和緯度相關(guān)性分析

青藏高原地區(qū)地形起伏大、氣候變化多端,極大地影響了Tm的準(zhǔn)確計算．文獻[17]表明Tm與Ts、高度和緯度具有較大的相關(guān)性．雖然Tm與Ts的線性關(guān)系已被證實,但在一些地區(qū)有出現(xiàn)Tm與Ts的關(guān)系不密切的現(xiàn)象[18],因此有必要分析在青藏高原地區(qū)Tm與Ts的關(guān)系．為了更好地分析Tm隨Ts、高度和緯度的變化關(guān)系,本文選取了2017年青藏高原地區(qū)13個探空站數(shù)據(jù),利用數(shù)值積分法計算出每個站點的Tm值并對其與相應(yīng)測站的Ts、高度及緯度的相關(guān)性分別進行探究,結(jié)果如圖2所示．

圖2表明青藏高原地區(qū)的13個探空站點的Tm值與Ts、高度和緯度均存在近似的線性變化關(guān)系．其中Tm和Ts、Tm和緯度為正線性相關(guān),Tm和高度為負(fù)相關(guān)．Tm和Ts的相關(guān)性最強,相關(guān)系數(shù)為0.879,Tm和高度、緯度則次之．由此可見,考慮Tm隨Ts、高度及緯度的變化并對其進行相應(yīng)改正是準(zhǔn)確計算Tm的關(guān)鍵,在建立Tm模型時可以考慮加入這三個要素的變化．

2.2 TPTm模型的建立

文獻[5]已經(jīng)提出了一個簡單的線性關(guān)系式Tm=a+bTs來表示Tm和Ts之間的關(guān)系．然而,由于Tm同時也與高度和緯度有顯著的相關(guān)性,因此相比于簡單的線性關(guān)系式,建立一個顧及多種因素的多因子Tm經(jīng)驗?zāi)Ｐ涂梢愿_地獲取所需的Tm值．除空間因素外,Tm和Ts存在著時間變化,其表現(xiàn)為年周期和半年周期變化[17]．因此,在用Tm和Ts的簡單線性關(guān)系式計算Tm時,會引入周期殘差[19]．為了探討Tm殘差的季節(jié)變化,以2013—2018年玉樹探空站獲取的Tm和Ts為例,先建立其線性關(guān)系式,然后將實測的Tm與關(guān)系式得到的Tm相減求出其Tm殘差,結(jié)果如圖3所示．

由圖3可知,Tm殘差具有明顯的季節(jié)變化,主要體現(xiàn)為年周期和半年周期的變化特性．一個好的Tm模型應(yīng)該同時考慮空間和時間的變化,因此在構(gòu)建Tm模型時還需顧及Tm殘差的年周期和半年周期變化．基于上述內(nèi)容,本文對Tm和Ts線性模型進行修正,增加了調(diào)和函數(shù),并在模型表達式中考慮高度和緯度的變化．模型公式如下:

Tm(Ts,h,φ,DOY)=A0+A1Ts+A2h+A3φ+

(8)

式中,h為海拔高(km),φ為大地緯度,DOY為年積日,A0,A1,A2,A3,B1,B2,B3,B4均為模型系數(shù)．本文采用均勻分布在青藏高原地區(qū)的2014—2017年13個探空站數(shù)據(jù),利用最小二乘法對式(8)中的系數(shù)進行求解,由此得到青藏高原地區(qū)的顧及時空因素的大氣加權(quán)平均溫度新模型(TPTm模型),求解的模型系數(shù)如表2所示．

TPTm模型是一個適用于青藏高原地區(qū)的多因子局域經(jīng)驗?zāi)Ｐ?使用者僅需提供目標(biāo)處的Ts、海拔高、緯度和年積日,即可通過模型計算出目標(biāo)處的Tm值．此外,采用與TPTm模型同樣的建模數(shù)據(jù),將式(7)中的a,b系數(shù)進行重新計算,構(gòu)建一個青藏高原地區(qū)Bevis精化模型(Bevis-TP模型),求取的a,b值分別為83.2和0.66,后續(xù)將該模型與TPTm模型進行對比,以衡量TPTm模型中考慮時空因素所帶來的改進．

圖3 玉樹探空站Tm殘差時序變化Fig.3 Variations of the Tm residual errors from 2013 to 2018 in Yushu radiosonde station

表2 TPTm模型的系數(shù)

2.3 TPTm模型的精度驗證

本文采用偏差(bias)和均方根誤差(RMS)兩種指標(biāo)來作為模型的精度評估標(biāo)準(zhǔn),其表達式為

(9)

(10)

圖4 使用2018年探空站數(shù)據(jù)檢驗不同模型的偏差分布Fig.4 Distribution of biases for different models tested using radiosonde profiles of 2018

為了評估新建立的TPTm模型的精度,本文以2018年的青藏高原地區(qū)相應(yīng)的13個探空站數(shù)據(jù)作為參考值,來驗證TPTm模型的精度,同時與常用的Bevis模型、Bevis-TP模型和目前表現(xiàn)不錯的GPT2w對流層模型進行對比分析．TPTm模型和Bevis模型所需要的Ts值均由相應(yīng)的探空站提供．不同模型計算出的Tm值分別與參考值進行差值對比分析,由此可計算出2018年不同模型的年均偏差和RMS值,結(jié)果分別見表3、圖4和圖5．

由表3可知:Bevis模型表現(xiàn)出最大的正偏差值和年均偏差值,分別為6.45 K和4.68 K,GPT2w-5表現(xiàn)出顯著的負(fù)偏差值,其大小變化范圍從-12.1K至2.5 K,而TPTm模型具有最小的年均偏差值,為0.07 K,其偏差變化范圍也是4個模型里最小的,大小變化范圍從-1.35 K至1.07 K;對于RMS值,Bevis表現(xiàn)出最大的RMS值,為6.07 K,而TPTm模型的年均RMS值最小,為2.76 K．局域精化后的Bevis-TP模型比Bevis模型在精度上有明顯的提升,表明基于本地數(shù)據(jù)構(gòu)建的Bevis-TP模型相比外地數(shù)據(jù)構(gòu)建的Bevis模型可靠性更高．相比于GPT2w-5模型,GPT2w-1模型的偏差值和RMS值都有明顯改善,沒有出現(xiàn)過大的絕對偏差值和RMS值,表明水平分辨率較高的GPT2w-1模型相比水平分辨率低的GPT2w-5其精度有一定改善且相對穩(wěn)定．相比于Bevis、Bevis-TP、GPT2w-5和GPT2w-1模型,TPTm模型的精度分別提高54.5%(3.31 K)、30.8%(1.23 K)、36.3%(1.57 K)和27.6%(1.05 K),表明TPTm模型在青藏高原地區(qū)計算的Tm值明顯比其他模型的精度高,特別是對Bevis模型其精度有較好的改善．

表3 使用2018年探空站數(shù)據(jù)驗證不同模型的精度對比

圖5 使用2018年探空站數(shù)據(jù)檢驗不同模型的RMS誤差分布Fig.5 Distribution of RMS errors for different models tested using radiosonde profiles of 2018

由圖4可以看出,Bevis模型在青藏高原地區(qū)多個站點存在顯著的正偏差,GPT2w-5和GPT2w-1模型則存在明顯的負(fù)偏差,而TPTm模型在青藏高原地區(qū)整體表現(xiàn)出較小的偏差和良好的穩(wěn)定性．由圖5可知,Bevis模型在青藏高原地區(qū)的RMS誤差整體偏大．其中在青藏高原地區(qū)的東部,Bevis模型的RMS值明顯大于其他地區(qū),且這些誤差大的站點位于海拔較高的地區(qū),表明在對于像青藏高原那樣的高海拔地區(qū),Bevis模型不能準(zhǔn)確地計算當(dāng)?shù)氐腡m值．而利用青藏高原局地的探空資料構(gòu)建的Bevis-TP模型精度表現(xiàn)相比Bevis模型有明顯的改善,尤其在青藏高原東部高原地區(qū)改善明顯．這是由于Bevis模型是利用美國本土(27°～65°N地區(qū))探空資料建立的,該模型在青藏高原地區(qū)使用時可能會產(chǎn)生局域模型系統(tǒng)誤差,使獲得的Tm值并非最優(yōu)．因此,在有氣象資料條件下最好采用當(dāng)?shù)貧庀髷?shù)據(jù)來構(gòu)建Tm模型．但Bevis-TP模型與TPTm模型相比,其偏差和RMS誤差相對較大．除了青藏高原地形起伏較大的原因外,Bevis和Bevis-TP模型沒有考慮Tm和Ts的季節(jié)變化也是導(dǎo)致其誤差大的主要因素之一．GPT2w模型與Bevis模型相比其表現(xiàn)相對較好,但還存在個別偏大的RMS誤差．GPT2w模型雖然考慮了年周期和半年周期變化,但由于忽略了Tm隨高程的變化[20],從而導(dǎo)致其在地形起伏較大的青藏高原地區(qū)使用時出現(xiàn)顯著的系統(tǒng)誤差．且GPT2w模型提供逐日的Tm值,無法滿足獲取逐小時Tm值的需求．然而,TPTm模型克服了這些缺點,它在建模時同時考慮了Ts、地理位置(海拔和緯度)和時間的變化,只要用戶提供逐小時的Ts值,就可得到對應(yīng)每小時的Tm值．與其他模型相比,TPTm模型的總體RMS誤差較小,且相對穩(wěn)定,特別是在高海拔地區(qū)表現(xiàn)良好．這些評估結(jié)果表明TPTm模型在青藏高原地區(qū)具有較好的性能,表現(xiàn)出較高的精度和良好的穩(wěn)定性．

2.4 Tm對GNSS-PWV計算的影響

建立青藏高原地區(qū)的新Tm模型是為了提高Tm的計算精度,其最終目的是提高GNSS反演PWV的精度．但是一般情況下,GNSS基準(zhǔn)站與探空站不在同一地址,且大多數(shù)GNSS基準(zhǔn)站主要用于大地測量研究,未安裝氣象傳感器,所以難以全面、可靠地研究Tm對GNSS-PWV計算的影響．為此,本文采用文獻[17]提出的計算Tm對GNSS-PWV的影響的方法,并對計算結(jié)果進行分析:

(11)

圖6 利用2018年探空站資料檢驗Bevis模型、Bevis-TP模型、GPT2w模型及TPTm模型計算PWV的理論RMS誤差和相對誤差Fig.6 Theoretical ERMS,PWV and ERMS,PWV/VPW errors from the Bevis,Bevis-TP,GPT2w and TPTmmodels tested using radiosonde data of 2018

式中:ERMS,PWV是PWV的RMS誤差;ERMS,Π為轉(zhuǎn)換系數(shù)Π的RMS誤差;ERMS,Tm是Tm的RMS誤差;ERMS,PWV/VPW定義為PWV的相對誤差．Tm和PWV采用年均值計算．各模型在13個探空站上的ERMS,PWV和ERMS,PWV/VPW的計算結(jié)果如圖6所示．

從圖6可知,Bevis模型和GPT2w-5模型在多數(shù)站點上均表現(xiàn)出較大的ERMS,PWV和ERMS,PWV/VPW值,ERMS,PWV最大值分別為0.28 mm和0.56 mm,ERMS,PWV/VPW最大值分別為2.67%和4.86%．而TPTm模型的整體ERMS,PWV值均小于0.15 mm,平均ERMS,PWV值為0.10 mm．TPTm模型的ERMS,PWV/VPW平均值為1.02%,范圍為0.67%至1.27%,比其他模型更穩(wěn)定且更?。虼?TPTm模型計算出來的Tm值所帶來的誤差對計算出的GNSS-PWV的影響,相比其他模型要小且較穩(wěn)定,可以為青藏高原地區(qū)的GNSS反演出精確的PWV值提供更為準(zhǔn)確的Tm信息．

3 結(jié)語

1)本文采用2014—2017年13個探空站數(shù)據(jù),建立了顧及Ts、高度、緯度和季節(jié)變化的青藏高原地區(qū)大氣加權(quán)平均溫度模型(TPTm模型),并以2018年探空站獲取的Tm作為參考值,模型的平均偏差和平均RMS分別為0.07 K和2.76 K．

2)通過將TPTm模型與常用的Bevis模型、局域精化后的Bevis模型和GPT2w模型相比較可知,TPTm模型相比其他模型在青藏高原地區(qū)的平均偏差和RMS更小、精度更高,對青藏高原地區(qū)有著更好的適應(yīng)性、穩(wěn)定性．

3)從理論上分析了各模型計算的Tm對GNSS-PWV估計的影響,得到TPTm模型的ERMS,PWV和ERMS,PWV/VPW的平均值分別為0.10 mm和1.02 %,比其他模型更小且更穩(wěn)定．

總之,TPTm模型得到的Tm精度可靠,可以滿足青藏高原地區(qū)GNSS水汽探測的要求,從而進一步改善青藏高原地區(qū)GNSS大氣水汽反演的精度,更好地服務(wù)于青藏高原地區(qū)氣象學(xué)的研究工作．