連達軍,何琦敏,2,李 黎,2,富爾江,張克非
(1.蘇州科技大學 地理科學與測繪工程學院,江蘇 蘇州 215009;2.蘇州科技大學 北斗導航與環境感知研究中心,江蘇 蘇州 215009;3.北星空間信息技術研究院(南京)有限公司,江蘇 南京 211800;4.中國礦業大學 環境與測繪學院,江蘇 徐州 221116)
天頂對流層延遲(Zenith Tropospheric Delay,ZTD)和大氣可降水量(Precipitable Water Vapor,PWV)是改善氣象預報的重要大氣參數,同時也是影響全球導航衛星系統(Global Navigation Satellite System,GNSS)高精度定位的主要誤差源。ERA5是歐洲中期天氣預報中心提供的第五代大氣再分析數據集,能夠提供全球范圍內的小時氣象數據。由于ERA5數據具有較高的時空分辨率和精度,常用于構建ZTD等大氣參數的經驗模型[1-2],并作為GNSS和衛星遙感等水汽觀測技術的檢校源[3]。由于ZTD和PWV作為輸入數據源參與氣象預報時應滿足要求相應的精度閾值要求[4],國內外已有部分學者對ERA5資料反演的ZTD(ERA5-ZTD)和PWV(ERA5-PWV)進行了一些初步的評估。在全球范圍內,YU等[5]利用1 300個GNSS監測站反演的ZTD(GNSS-ZTD)對ERA5-ZTD精度進行了評估,兩者的整體偏差約為16.9 mm,全球范圍內精度不均勻,與站點的高程具有相關關系;上官明等[6]基于分析了基于全球180個GNSS站和180個探空站(Radiosonde Station,RS)反演的ZTD對ERA5-ZTD的精度進行了評估,其結果具有顯著的季節性特征,在夏季時的精度最低;ZHANG等[7]檢驗了ERA5、ERA- Interim(ERA5的上一代產品)、GNSS和RS反演的PWV數據集的一致性,以GNSS-PWV和RS-PWV作為參考值,ERA5-PWV的均方根誤差(Root Mean Square Error,RMSE)分別為1.8 mm和2.7 mm。此外,有不少學者分析了ERA5產品在我國省市區域范圍的適用性,相較于其他高精度大氣探測方法獲得的ZTD和PWV產品,ERA5-ZTD和ERA5-PWV的平均偏差分別在[15 mm,20 mm]和[1.5 mm,3 mm]區間范圍[8-12]。但上述的研究都是將兩者的評估分開,且時間尺度較短,因此得到的結果未能客觀反映ZTD和PWV兩者的整體精度關系。
此外,由于在極端天氣事件演變過程中,ZTD和PWV等大氣參數通常出現顯著的異常變化,因此ZTD和PWV可直接應用于部分區域性極端天氣的預警預報[13]。目前ZTD和PWV在氣象預報方面的應用主要包括降雨、臺風和對流風暴等方面[14],已能夠實現定性和部分定量的極端天氣短臨預警。本文利用中國區域2017-2019年間的14個GNSS站和89個RS站點的GNSS-ZTD和RS-PWV分別對ERA5-ZTD和ERA5-PWV進行精度評估。以2019年我國遭受的超強臺風“利奇馬”事件為例,分析“利奇馬”期間的ZTD和PWV時空分布及運移規律,構建臺風水汽移動因子,研究水汽和臺風的響應關系,為我國預報臺風等災害性天氣提供新的參考思路。
ZTD可根據大氣折射率在天頂方向的路徑積分表達式計算:

(1)
式中:h0和htro分別為地面和對流層頂高度(m),N為大氣折射率,dh為天頂方向路徑微分(m)。N難以用顯式的表達式計算,可用分層的延遲量近似表示[15]:
(2)
式中:Ni和Δhi分別為分層的平均大氣折射率和高度差(m),Pi、ei、和Ti分別為分層的平均總氣壓(hPa)、水汽分壓(hPa)和溫度(K),n為總分層數,k1,k2和k3均為常數,其值分別為77.6、64.8和3.776×105。
在式(2)中,參數T和P可直接從ERA5數據集獲得,e需要通過比濕和氣壓間接計算得到[16]:
(3)
式中:Q為比濕,可從ERA5數據集中直接獲得。
PWV可根據比濕在天頂方向的氣壓積分表達式得到,采用分層氣壓的方式計算[17]:
(4)
式中:Qi和ΔPi分別分層的平均比濕和氣壓差,g為當地重力加速度(m/s2)。
此外,由于ERA5數據反演的ZTD和PWV數據是基于ERA5格網點位置,因此使用ERA5反演的數據產品與GNSS和RS的結果對比時,還需考慮ERA5格網點與GNSS和RS的站點空間差異。其中,在水平方向上采用距離站點最近的四個格網點進行反距離加權的方法,獲得站點所在水平位置的ERA5-ZTD/PWV,然后在垂直方向上采用樣條曲線方法插值得到站點所在高度的ERA5結果。此外,由于ERA5和GNSS存在高程系統上的差異,其中ERA5數據產品是基于位勢高系統,而GNSS數據產品基于大地高系統,通過文獻[18]中的方法可實現兩者轉換。
我國東南沿海地區常年受到臺風的侵襲,是世界上臺風活動最活躍的地區之一。為了綜合評估ERA5-ZTD/PWV產品在我國區域的可用性,并分析ZTD和PWV的時空變化與臺風天氣事件之間的響應機理。本實驗所需的數據集主要包括如下:
1)2017-2019年ERA5位勢、溫度和比濕等參數的氣壓分層格網數據,空間分辨率、時間分辨率和垂直分辨率分別為0.25°×0.25°、1 h和37層(https://cds.climate.copernicus.eu/#!/ home)。
2)國際GNSS服務(International GNSS Service,IGS)提供的2017-2019年中國區域14個GNSS站點事后處理的ZTD產品,時間分辨率為5 min(https://cddis.nasa.gov/archive/gnss/products/troposphere/zpd/)。
3)懷俄明大學大氣科學系提供的2017-2019年中國區域89個探空站數據,時間分辨率為12 h(http://weather.uwyo.edu/upperair/sounding.html)。
4)中央氣象臺提供的2019年臺風最佳軌跡數據集[19-20],時間分辨率為3 h或6 h (http://tcdata.typhoon.org.cn)。
數據集(1)主要用于反演ERA5-ZTD和ERA5-PWV,數據集(2)和數據集(3)分別用于評估ERA5-ZTD和ERA5-PWV的精度,數據集(4)主要用于提取2019年“利奇馬”臺風的移動路徑,并計算路徑上的移動速度,分析ZTD和PWV時空變化與臺風臨近前后的響應關系。其中,數據集(2)中的14個IGS站點的大地坐標和高程見表1,數據集(2)的IGS站點和數據集(3)的探空站點分布見圖1。

圖1 實驗中GNSS站和探空站地理分布

表1 GNSS站點大地坐標及高程
IGS采用事后的精密衛星軌道和鐘差處理策略,能夠提供時間分辨率為5 min,且精度優于5 mm的ZTD產品[21]。其中在2017-2019年時間段內,中國區域范圍共有14個IGS站點的ZTD數據產品可用。因此,為了檢驗ERA5-ZTD的可用性,利用ERA5格網數據首先獲得所在的IGS站點ZTD值,以IGS提供GNSS-ZTD產品作為參考值,計算ERA5-ZTD在各個站點上的平均偏差(Bias)、RMSE和標準差(Standard Deviation,STD),統計結果見表2。

表2 14個GNSS站點的ERA5-ZTD統計結果
由表2結果可知,14個站點的ERA5-ZTD平均Bias、RMSE、STD分別為-6.3 mm、19.3 mm和17.5 mm,絕大部分測站精度在13.9~21.1 mm之間;其中誤差最大的站點為LHAZ,平均Bias、RMSE、STD分別為7.6 mm、34.3 mm和33.4 mm。分析LHAZ站點的ERA5-ZTD誤差偏大的主要原因可能是由于該站點位于氣象站點和數據匱乏的青藏地區,而ERA5的大氣數據產品是基于多種氣象數據的同化模型結果。因此在原始數據缺失的地區,同化后的氣象數據結果精度較差,而高海拔的青藏地區環境和條件使搭建氣象站變得困難,且難以長期維護,因此基于ERA5數據反演ZTD產品整體精度也將偏低。
探空氣球是一種高精度的高空大氣探測技術,能夠獲取探空站上空的多種氣象參數(例如:溫度、氣壓和相對濕度等)廓線信息,基于該廓線信息可以反演得到PWV。由探空站反演的PWV精度通常優于1.5 mm,但時間分辨率較低,通常為1 d 2次或3次。因此,RS-PWV常用來檢校GNSS、衛星遙感和數值氣象預報模式等水汽探測方法的精度。
為了進一步分析ERA5-PWV在中國區域內的精度,利用2017-2019年間中國區域內89個探空站反演的RS-PWV數據作為參考值檢驗ERA5-PWV的精度,并計算了兩者的相關系數值(圖2)。在圖2中,由NW至SE方向精度明顯下降,整體呈現西部內陸精度高,東南沿海精度低的趨勢,ERA5-PWV的整體平均Bias、RMSE和STD分別為0.3 mm,2.9 mm和2.7 mm,最大值分別為-2.29 mm、6.65 mm和6.25 mm。ERA5-PWV和RS-PWV具有顯著的相關關系,超過97%站點的相關系數大于0.9,平均和最小相關系數分別為0.97和0.84。由ERA5-ZTD和ERA5-PWV在中國區域內的精度評估結果表明,兩者具有較高的精度,能夠進一步用于分析與極端天氣事件的響應關系。

圖2 89個RS站點的ERA5-PWV統計結果(Bias、RMSE、STD和相關系數)分布
為了進一步分析ZTD和PWV與臺風事件之間的響應關系,本實驗以2019年的“利奇馬”臺風事件為例,提取了2019年8月8日0時-18時“利奇馬”臺風的活動區域范圍的ERA5-ZTD和ERA5-PWV,并進行了空間插值,獲得了PWV和ZTD的空間分布圖(圖3和圖4)。由圖可知,“利奇馬”臺風期間的ZTD和PWV空間分布非常相似。當臺風內部攜帶了大量水汽,臺風路徑上有一個近似的旋轉黃色圓(高PWV)能夠顯著表征當時臺風覆蓋區域和所處位置。而ZTD是濕延遲和干延遲的綜合,結果中,臺風期間的ZTD顯著增大,可能是由于水汽的增加引起的濕延遲增大。此外,從圖中PWV和ZTD移動趨勢上可知,PWV和ZTD的移動方向具有明顯東南向西北方向趨勢,與“利奇馬”臺風的路徑方向是保持一致的。

圖3 2019年8月8日0時-18時每6 h的PWV空間分布與“利奇馬”臺風路徑

圖4 2019年8月8日0時-18時每6 h的ZTD空間分布與“利奇馬”臺風路徑
因此,為了進一步研究“利奇馬”臺風期間的水汽時序特征,在路徑周圍選取了6個位于臺風中心附近ERA5格網點,分析臺風來臨前后的PWV變化。同時,為了對比多個時間序列的特征變化,對原始PWV序列進行標準化:
(5)


圖5 2019年8月6日0時-8月9日0時“利奇馬”臺風路徑周圍所選的ERA5格網點與1 h SPWV時間序列
由圖5可知,由于臺風攜帶了大量水汽,因此在來臨過程中,SPWV會出現一段時間的增長,之后達到最大值,隨著臺風的逐漸遠離,SPWV又開始呈現下降的趨勢。同時,由于臺風均經過上述6個ERA5點,因此圖5b中的6個SPWV序列非常相似,均能夠找到臺風遠離時SPWV的最大值點。但由于臺風來臨的先后次序不同,因此在上述幾個SPWV序列中最大值點出現的時間存在偏移,且該時間偏移方向和臺風到達的路徑方向也是一致的。由圖3可知,水汽的移動方向和臺風移動方向是接近的,為了進一步定量分析水汽的移動速度和臺風速度關系,假設臺風內的水汽首先到達了第1個格網點,且正在往附近第2個格網點方向移動,示意圖見圖6。

圖6 臺風水汽的平面運動模型
圖6中,d12為沿著臺風路徑方向行走的路程,l12為兩格網點之間的距離。因此,當臺風離開1號點到達2號點,所經過的時間為t12,即有:
d12≈l12cos(θ-θ12)。
(6)
式中:θ和θ12分別為水汽的移動方向和l12的方向角。不考慮水汽的收支和狀態變化等影響,當移動至最后一個格網點(圖5中6號點),即有:
(7)
式中:d為沿著臺風路徑方向行走的路程總長,di(i+1)、li(i+1)和θi(i+1)為分段的路程、格網點距離和方位角。
假設水汽的移動速度v保持不變,由圖5b可以得到各相鄰格網點的水汽到達時間差ti(i+1),由此可以通過下式計算臺風運動過程中的水汽平均移動速度。
(8)
根據上述公式即計算得到“利奇馬”臺風在此期間的水汽平均移動速度為12.4 km/h,而通過臺風的軌跡數據計算的臺風移動速度為16.1 km/h。兩者的計算結果非常接近,但由于臺風內的水汽還存在過程損失等情況,因此在數值上水汽的速度將小于實際的臺風移速,但仍可以作為一個重要參考因子對臺風的運動狀態進行預報。
為了研究ERA5數據反演的ZTD和PWV在中國區域內的可用性,本文利用2017-2019年中國區域內的14個IGS站和89個RS站估計的ZTD和PWV分別對ERA5-ZTD和ERA5-PWV進行精度評估。結果表明,中國區域內的ERA5-ZTD和ERA5-PWV精度較高,ERA5-ZTD的平均Bias、RMSE、STD分別為-6.3 mm、19.3 mm和17.5 mm,除高海拔的青藏區域LHAZ站存在較大偏差外,其余站點的精度均在13.9~21.1 mm范圍內;ERA5-PWV的平均Bias、RMSE和STD分別為0.3 mm,2.9 mm和2.7 mm,精度分布整體呈現西部內陸精度高和東南沿海低的趨勢。
此外,本文還分析ZTD和PWV與臺風事件的響應關系,以2019年“利奇馬”臺風為例,分析了臺風期間ZTD和PWV的空間分布和時間序列。結果表明,ZTD和PWV具有相似的分布規律,移動方向呈現與臺風的路徑方向一致趨勢,當臺風臨近時均出現顯著的增長,達到最大值隨臺風遠離開始下降。因此,為定量研究水汽的移動速度和臺風移動關系,提出了使用臺風路徑附近的ERA5格網點水汽序列計算水汽速度的方法,計算了“利奇馬”臺風在經過6個ERA5格網點的平均水汽移動速度為12.4 km/h,與臺風在該段路徑上的平均移動速度16.1 km/h非常接近。以上的研究表明,ERA5-ZTD和ERA5-PWV在中國區域范圍內具有較高的精度,ZTD和PWV可以作為臺風預警預報研究的重要輔助氣象因子。本文所使用的ERA5水汽產品是事后解算的結果,隨著實時GNSS等水汽探測技術發展,基于多源水汽的實時臺風監測技術將得到進一步發展。