李慧, 張亞豪, 朱丹彤, 高旭昂

(華北水利水電大學(xué),河南 鄭州 450046)

隨著空間大地測量技術(shù)的發(fā)展,基于地基全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System,GNSS)的空間探測技術(shù)自1992年起受到廣泛關(guān)注[1]。GNSS信號在傳播過程中受對流層影響出現(xiàn)路徑彎曲和速度衰減,使GNSS觀測值受到對流層延遲(Zenith Tropospheric Delay,ZTD)的影響。對GNSS觀測值進(jìn)行高精度解算得到ZTD估計值,結(jié)合多種氣象數(shù)據(jù)可反演出測站上空的大氣可降水含量(Precipitable Water Vapor,PWV)。傳統(tǒng)的PWV觀測方法有無線電探空儀、地基微波輻射計等。無線電探空儀由于其時空分辨率低,難以滿足高時空分辨率的PWV監(jiān)測需求。地基微波輻射計由于其空間分辨率低且價格昂貴無法被普及。相比于其他手段,地基GNSS水汽探測具備全天候、高精度、高時間分辨率和均一性好等優(yōu)勢,近些年已成為大氣水汽觀測的一種重要手段[2],并成功應(yīng)用到降雨預(yù)報[3]、水汽三維分布[4]、厄爾尼諾[5]、長期氣候變化[6]等研究中。如:單路路等[7]以PWV變化率、每月PWV閾值和PWV變化量3個因素作為參考要素,提出了一種利用GNSS-PWV進(jìn)行短期臨時降雨預(yù)測的方法;ZHAO Q Z等[8]基于月平均PWV、季節(jié)性PWV/ZTD的變化和它們的一階導(dǎo)數(shù)等5個預(yù)測因子,提出了一種改進(jìn)的降雨預(yù)報模型;ZHENG Y X等[9]結(jié)合機(jī)器學(xué)習(xí)構(gòu)建了一種新型的GNSS可降水量實時遙感預(yù)測模型。

河南省位于我國中原地區(qū),北、西、南三面環(huán)山,地勢西高東低。夏季,海洋上的水汽在夏季風(fēng)和副熱帶高壓的驅(qū)動下源源不斷向內(nèi)陸傳送,受地形阻攔,水汽在河南上空積累,并與河南當(dāng)?shù)氐膮^(qū)域性對流系統(tǒng)疊加,極易發(fā)生暴雨事件。2021年7月17—23日河南省遭遇歷史罕見特大暴雨,據(jù)國家自然災(zāi)害災(zāi)情管理系統(tǒng)統(tǒng)計,本次暴雨共造成河南省150個縣(市、區(qū))、1 616個鄉(xiāng)鎮(zhèn)、1 391.28萬人受災(zāi),同時引發(fā)了嚴(yán)重的城市內(nèi)澇、農(nóng)田積澇,造成了嚴(yán)重?fù)p失。大氣水汽是大氣中最為活躍的成分,是降雨發(fā)生的主要驅(qū)動因子之一。因此,分析大氣水汽在暴雨期間的時空演變特征對于進(jìn)一步理解暴雨發(fā)生機(jī)制、推動降雨預(yù)報現(xiàn)代化服務(wù)體系發(fā)展,具有十分重要的科學(xué)意義。鑒于以往有關(guān)鄭州“7·20”特大暴雨事件的研究都集中于降雨機(jī)制的分析,本文從大氣水汽的角度出發(fā),基于國際GNSS服務(wù)組織(International GNSS Service,IGS)分析中心的事后星歷、鐘差產(chǎn)品等,采用精密單點定位(Precise Point Positoning,PPP)解算GNSS觀測值,結(jié)合ERA-5氣象觀測資料得到時間分辨率為5 min的大氣可降水含量,對河南“7·20”特大暴雨期間大氣水汽與降雨量進(jìn)行對比分析,得到大氣水汽與降雨量之間的協(xié)同變化特征。

1 數(shù)據(jù)與方法

1.1 GNSS大氣水汽

1.1.1 GNSS大氣水汽反演基本流程

1)利用非差非組合PPP,結(jié)合IGS事后觀測星歷處理CORS網(wǎng)觀測數(shù)據(jù)。對于電離層延遲,將其作為白噪聲進(jìn)行參數(shù)估計,并把對流層延遲ZTD、整周模糊度、接收機(jī)和衛(wèi)星鐘差等視為未知參數(shù),通過多參數(shù)估計法建立的誤差方程和法方程進(jìn)行觀測值解算,求得對流層延遲ZTD。

2)結(jié)合ERA-5氣象觀測資料以及天頂靜力延遲計算模型得到天頂靜力學(xué)延遲(Zenith Hydrostatic Delay,ZHD),進(jìn)而得到天頂濕延遲(Zenith Wet Delay,ZWD)。

3)基于ERA-5氣象觀測資料計算大氣加權(quán)平均溫度,確定水汽轉(zhuǎn)換系數(shù),將天頂濕延遲ZWD轉(zhuǎn)換為可降水汽含量PWV[10]。

1.1.2 非差非組合模型

地基GNSS觀測值解算方法有雙差網(wǎng)解法和精密單點定位法兩種。對于雙差網(wǎng)解法,當(dāng)不同測站之間的接收機(jī)距離較近時,接收機(jī)上空的對流層延遲存在較強(qiáng)的相關(guān)性,僅能夠獲得接收機(jī)和對流層延遲之間的相對值,此時必須引入長距離測站消除這種強(qiáng)相關(guān)性才能獲得各個測站的對流層延遲;而精密單點定位方法則無需顧及GNSS測站之間的相關(guān)性,可以直接求得各接收機(jī)的對流層延遲,且相較于雙差網(wǎng)解法的計算效率更高。非差非組合模型作為PPP的一種解算方案,在處理多頻GNSS數(shù)據(jù)時更加靈活,可以避免無電離層組合對噪聲的放大效應(yīng)。因此,采取非差非組合PPP進(jìn)行GNSS測站觀測數(shù)據(jù)解算。非差非組合觀測模型采用原始偽距和載波相位觀測方程,將各衛(wèi)星斜路徑上的電離層延遲作為參數(shù)進(jìn)行估計。原始偽距和載波相位觀測量的計算式分別為:

Pi=ρ+c(Δtr-Δts)+Tr+aiI+mi+εpi;

(1)

Li=ρ+c(Δtr-Δts)+Tr-aiI+Bi+Mi+εLi。

(2)

其中,

式中:Pi、Li分別為偽距和載波相位觀測量;i為信號頻率,i=1、2;ρ為衛(wèi)星到接收機(jī)相位中心的幾何距離;c為光在真空中的傳播速度;Δtr、Δts分別為修正后的接收機(jī)鐘差和衛(wèi)星鐘差;Tr為對流層延遲;I為修正的f1頻率信號電離層延遲;mi和Mi分別為偽距和載波相位觀測量的多路徑誤差;Bi為修正的模糊度參數(shù);εpi和εLi分別為偽距和相位觀測噪聲;DCBr、DCBs分別為接收機(jī)端硬件碼延遲、衛(wèi)星端硬件碼延遲;γi、ai為頻率系數(shù)。

1.1.3 GNSS對流層延遲

GNSS衛(wèi)星信號傳播到地面接收機(jī)的過程中,由于穿過對流層發(fā)生折射導(dǎo)致信號路徑發(fā)生彎曲、信號傳播速度發(fā)生變化,進(jìn)而產(chǎn)生對流層延遲[11]。信號延遲量ΔL可由式(3)計算出。

(3)

(4)

式中:L為信號傳播距離;N為信號傳播路徑上的大氣折射率;T為氣溫,℃;Pd和Pv分別為干空氣壓強(qiáng)和水汽壓強(qiáng),hPa;Zd和Zv分別為干、濕壓縮因子;k1、k2、k3為氣體常數(shù),其值分別取77.604、64.290、3.776×105K2·hPa-1。

令N=Nh+Nw,Nh和Nw分別表示大氣折射的靜力學(xué)部分和非靜力學(xué)部分的折射率,將N代入公式(3)可得:

(5)

式中ΔLh、ΔLw分別為對流層的干、濕延遲量。信號傳播路徑與其對流層延遲量一一對應(yīng),然而在實際的GNSS數(shù)據(jù)處理中,將所有的GNSS信號路徑上的延遲量(ΔL)都作為未知數(shù)進(jìn)行解算會導(dǎo)致法方程的秩虧問題。因此,在實際解算中采用映射函數(shù)將信號傳播路徑上的對流層延遲量轉(zhuǎn)化為對流層天頂總延遲量(Zenith Tropospheric Delay,ZTD),并將其作為未知值進(jìn)行平差計算。ZTD的表達(dá)式[12]為:

ZTD=ZHD+ZWD。

(6)

式中:ZTD為由非偶極分量引起的天頂總延遲量;ZHD為由非偶極分量引起的天頂干延遲量;ZWD為由偶極分量引起的天頂濕延遲量。

ZHD的解算精度可達(dá)到亞毫米級,而天頂濕延遲量ZWD的最高精度僅達(dá)到厘米級,這顯然無法滿足氣象實際應(yīng)用的要求[12]。對于GNSS-PWV來說,通常是由GNSS高精度數(shù)據(jù)處理得到毫米級精度的天頂總延遲量ZTD,再利用Saastamoinen模型結(jié)合測站位置的氣象信息計算得到ZHD值,最后得到ZWD估計值?；赟aastamoinen模型的ZHD計算公式如下:

(7)

f(Lat,SH)=1-0.002 6cos(2Lat)-0.000 28SH。

(8)

式中:Ps為GNSS測站處的地面大氣壓,hPa;f(Lat,SH)是測站緯度Lat(單位為rad)和大地高SH(單位為km)的函數(shù)。

1.1.4 轉(zhuǎn)換系數(shù)的確定

ZWD和PWV之間的轉(zhuǎn)換關(guān)系[13]為:

PWV=Π·ZWD;

(9)

(10)

式中:Π為和對流層溫度有關(guān)的轉(zhuǎn)換系數(shù);ρw為液態(tài)水的密度,ρw=1 000 kg·m-3;Tm為加權(quán)平均溫度,K;Rv為水汽氣體常數(shù),Rv=461.51 J/(K·kg);k′2為氣體常數(shù),k′2=17±10 K/hPa。

加權(quán)平均溫度Tm是轉(zhuǎn)換系數(shù)計算中最重要的參數(shù)之一,本文采用張洛愷等[15]提出的適用于鄭州區(qū)域的Tm模型進(jìn)行計算。該模型是基于測站的地表溫度Ts、大氣壓強(qiáng)P和地面水汽壓Pv的多因子模型。Tm的計算式為:

Tm=0.467Ts-0.239P+0.018Pv+380.88。

(11)

本文采用的數(shù)據(jù)為河南地區(qū)CORS網(wǎng)中8個GNSS測站“7·20”暴雨發(fā)生前后2021年7月17—23日(由于CHSM測站缺少17日數(shù)據(jù),因此采用18—23日)共7 d 的觀測數(shù)據(jù)以及河南鄭州氣象站50783號19—23日的逐時降雨量數(shù)據(jù)進(jìn)行精度分析。其中GNSS觀測值的時間分辨率為5 min,GNSS-PWV反演過程中所涉及的氣象參數(shù)采用ERA-5再分析數(shù)據(jù)(小時分辨率)。為消除兩者之間的時間偏差,采用線性插值的方法獲取5 min分辨率的氣象數(shù)據(jù)。所選GNSS測站的空間分布、測站信息和氣象站站點信息分別如圖1所示并見表1。

表1 河南地區(qū)GNSS測站和氣象站信息

圖1 河南省GNSS測站分布

1.2 ERA-5大氣水汽

1.2.1 ERA-5大氣水汽反演原理

ERA-5是由歐洲中期天氣預(yù)報中心發(fā)布的第五代全球氣象再分析資料,可以提供從1979年到現(xiàn)在的氣象再分析數(shù)據(jù)。ERA-5分層及地表數(shù)據(jù)的時間分辨率為1 h,水平方向的空間分辨率為31 km(0.25 °),高程方向?qū)⒌孛嬷?0 km高空的空間分成37個氣壓層(從1 000 hPa至1 hPa),分層數(shù)據(jù)包含溫度、氣壓、比濕、風(fēng)速、位勢等數(shù)據(jù),地表數(shù)據(jù)有氣壓、2 m露點溫度、位勢及比濕[16]。

基于ERA-5再分析數(shù)據(jù),通過對地面上空每個氣壓層中的氣象參數(shù)進(jìn)行數(shù)值積分可獲取目標(biāo)點上空的PWV,具體計算公式為:

(12)

其中:K為大氣廊線上的水壓層的個數(shù);g為重力加速度,m/s2;q為比濕,kg/kg;P為壓強(qiáng),hPa。

由于無法獲得大氣廓線的連續(xù)數(shù)據(jù),故采用各個等壓面上氣象參數(shù)的差分值代替積分進(jìn)行PWV計算[16],即:

(13)

ΔPi=Pi+1-Pi。

式中:ΔPi為相鄰兩個等壓面的差值;qi為相鄰等壓面的平均比濕;Pi+1、Pi分別為第i+1、i個等壓面的壓強(qiáng)。對于比濕q,先采用各個等壓面的溫度數(shù)據(jù)計算等壓面上的飽和水汽壓,再利用飽和水汽壓與比濕的關(guān)系計算比濕數(shù)據(jù)。計算公式分別為:

(14)

(15)

(16)

式中:Psa為飽和水汽壓;T為開氏溫度;RH為相對濕度。

1.2.2 高程改化

采用再分析資料中大氣廓線上的相對濕度、溫度和壓強(qiáng)依次計算各個等壓面位置處的比濕,采用式(13)求得任意等壓面上空的PWV數(shù)據(jù),但是由于GNSS測站的高程與等壓面上的值并不相同,并且各個高程系統(tǒng)并不統(tǒng)一。因此,為了計算GNSS測站上空的PWV數(shù)據(jù),需要先將不同高程基準(zhǔn)下的參數(shù)統(tǒng)一到同一高程基準(zhǔn)下,然后將等壓面上的氣象參數(shù)進(jìn)行高度改化,將其歸算到指定高度處。

GNSS采用的是大地高系統(tǒng),再分析數(shù)據(jù)采用的是位勢系統(tǒng),本文以正高系統(tǒng)為基準(zhǔn)。對于再分析數(shù)據(jù),首先,利用位勢計算各個等壓面的位勢高。其計算公式為:

(17)

式中:GPH為位勢高;GP為位勢;gLat為等壓面的重力加速度。為了便于計算,本文忽略重力加速度的垂直變化,直接采用下式計算gLat。

gLat=9.806 2[1-2.644 2×10-3cos(2Lat)+
5.8×10-6cos2(2Lat)]。

(18)

然后,將再分析數(shù)據(jù)的位勢高轉(zhuǎn)化為正高。其計算公式為:

(19)

(20)

式中:H為正高;RLat為地球有效半徑。

在氣象參數(shù)統(tǒng)一轉(zhuǎn)化到正高系統(tǒng)下后,為了獲取指定高度處的氣象參數(shù),需要進(jìn)行高程改正以消除高程非一致性導(dǎo)致的偏差。若目標(biāo)點高于大氣廓線的最底層,則基于上下兩層等壓面的氣象參數(shù),通過線性插值進(jìn)行高程改化;若目標(biāo)點的高程低于大氣廓線的最底層,則采用外推法進(jìn)行高程改化。

2 結(jié)果與分析

2.1 GNSS-PWV精度分析

本文利用ERA-5再分析資料反演得到的PWV(簡稱為ERA-PWV)對GNSS-PWV進(jìn)行精度評估。必須指出的是,由于ERA-5的時空分辨率與GNSS的不同,為此采用以下措施進(jìn)行處理:在空間域上,首先求得GNSS站點周圍4個ERA-5格網(wǎng)點位置上的PWV值,然后采用反距離加權(quán)法求得GNSS站點上空的PWV值;在時間域上,采用窗口長度為1 h的時間窗口法求取ERA-5歷元前后半小時的GNSS-PWV平均值,并以此平均值作為該歷元的GNSS-PWV并與ERA-PWV進(jìn)行比較。

精度比較采用的統(tǒng)計指標(biāo)分別為偏差(Bias)、均方根誤差(Root Mean Square,RMS)和相關(guān)系數(shù)(ρ)。偏差指的是期望輸出與真實數(shù)據(jù)的差別,用來描述模型對本訓(xùn)練集的擬合程度;均方根誤差衡量的是觀測量與真值的偏差;相關(guān)系數(shù)是衡量向量相似度的值,無量綱。偏差、均方根誤差、相關(guān)系數(shù)的計算公式分別為:

(21)

(22)

(23)

表2 河南各測站ERA-PWV和GNSS-PWV的偏差、均方根誤差和相關(guān)系數(shù)

圖3 河南地區(qū)各測站GNSS-PWV精度分布圖

由圖2可以看出:①8個測站的PWV值自7月18日起逐漸增大,并在7月20日達(dá)到峰值,此時GNSS-PWV與ERA-PWV的結(jié)果更加接近,二者的差值在±5 mm內(nèi)。②兩種PWV數(shù)據(jù)在55～60 mm范圍內(nèi)的偏差較大,其原因可能為18日暴雨發(fā)生前或者23日暴雨結(jié)束后PWV數(shù)據(jù)受到空氣流動、溫度、濕度等因素的影響;兩種PWV數(shù)據(jù)在60～70 mm范圍內(nèi)的離散程度較小,分布比較集中,說明兩者差值較小,此情況主要發(fā)生在20日降雨集中時,PWV數(shù)據(jù)受其他因素的影響較小。

圖3(a)反映了GNSS-PWV與ERA-PWV的均方根誤差的空間分布,可以明顯的看出,隨著緯度的增加其均方根誤差也明顯增加;圖3(b)反映了偏差、均方根誤差、相關(guān)系數(shù)的緯度分布,由此可以看出,偏差、均方根誤差、相關(guān)系數(shù)與經(jīng)度的增大并沒有明顯的關(guān)系。

由表2可知:偏差的最小值為-9.095 mm,最大值為6.483 mm,平均值為-0.731 4 mm,考慮到當(dāng)時暴雨天氣對偶然誤差的影響,總體偏差較小;均方根誤差最大值、最小值、平均值分別為2.689、1.795、2.290 mm。這表明GNSS-PWV整體具有較好的可信度。

圖4為河南地區(qū)8個測站一周內(nèi)GNSS-PWV與ERA-PWV的散點圖。

圖4 GNSS-PWV與ERA-PWV的散點圖

圖4表明GNSS-PWV與ERA-PWV具有很好的相關(guān)性,GNSS-PWV的時間分辨率較高,獲得的數(shù)據(jù)相對較多,具有較高的可靠性。

2.2 水汽分布特征分析

為了形象直觀地展現(xiàn)暴雨期間河南各測站上空PWV的變化過程,采用面積圖表示每個測站的GNSS-PWV與ERA-PWV,結(jié)果分別如圖5和圖6所示。由圖5和圖6可知:①7月17—18日PWV逐漸增加,直至7月20日和7月21日,PWV逐漸達(dá)到最大值,同時在7月21日之后PWV明顯下降。②7月21日之后雖然ZZXZ上空的PWV開始減小,但是XXYY、XXCY、JZJZ的還有不同程度的增加,但相比于7月20日的增幅較小,仍居于60 mm之上,說明7月20日之后的降雨態(tài)勢沒有減小,反而出現(xiàn)了持續(xù)的小高峰。③7月23日之后,PWV持續(xù)減少,說明降雨態(tài)勢逐漸減小且趨于穩(wěn)定。④整個一周的時間中,水汽一直處于很高的狀態(tài),這是因為受到臺風(fēng)“煙花”和“查帕卡”的影響,水汽不斷在測站上空聚集[20]。

圖5 河南各測站GNSS-PWV變化過程

圖6 河南各測站ERA-PWV變化過程

為進(jìn)一步分析降雨過程中測站所在地區(qū)降雨量與PWV變化的關(guān)系,采用折線圖表示GNSS-PWV(圖7)。圖7中的柱狀圖表示鄭州市氣象站實際降雨量信息。

圖7 河南各測站可降水含量(GNSS-PWV)變化過程

由圖7可知:①GNSS-PWV的變化趨勢與鄭州市氣象站上空逐小時降雨量均呈現(xiàn)先增加后降低的趨勢,并且水汽的峰值提前于降雨量的峰值,并在暴雨前3 h左右達(dá)到最大值,這與之前的研究結(jié)論相符合[26]。②從空間維度上看,河南省內(nèi)的GNSS-PWV由中部向四周逐漸降低,鄭州市GNSS-PWV的峰值最高,達(dá)到了70 mm/5 min,并且當(dāng)時鄭州中心及周邊區(qū)域出現(xiàn)了特大暴雨,降雨量全省最高,達(dá)到了200 mm/h。③從時間維度上看,河南8個測站的GNSS-PWV變化均呈現(xiàn)從低到高的狀態(tài),可以很好地解釋暴雨的變化過程。綜合以上分析可知,GNSS-PWV不但具有較高的時間分辨率,而且與實際降雨量有較高的符合程度,可以為極端暴雨短時間預(yù)報提供數(shù)據(jù)支撐。

3 結(jié)論

水汽作為大氣中重要的氣體成分之一,對大氣變化起著重要的影響。目前常規(guī)探測大氣水汽的方法有無線電探空技術(shù)、氣象觀測衛(wèi)星、氣象飛機(jī)探測等技術(shù),但常規(guī)技術(shù)存在成本高、精度低、探測范圍小等不足。相比于常規(guī)方法,GNSS數(shù)據(jù)具有高精度、高時間分辨率等特點,利用GNSS數(shù)據(jù)進(jìn)行水汽研究對于提供災(zāi)害性天氣的預(yù)報工作具有重要意義。本文基于河南省8個GNSS測站的觀測值研究了河南“7·20”特大暴雨期間的水汽分布特征。主要結(jié)論如下:

1)以河南省CORS網(wǎng)中8個GNSS測站在特大暴雨期間(7月17—23日)的觀測數(shù)據(jù)為例,解算得到GNSS-PWV數(shù)據(jù),并基于ERA-PWV評估了GNSS-PWV的精度。通過對比分析得到河南省大氣可降水含量在空間分布上隨著緯度的增加偏差逐漸增大。此外,根據(jù)降雨量統(tǒng)計數(shù)據(jù)發(fā)現(xiàn)GNSS-PWV的變化趨勢與實際降雨過程基本符合。

2)將GNSS-PWV與ERA-PWV進(jìn)行對比分析,發(fā)現(xiàn)其偏差主要分布在5 mm范圍之內(nèi),利用GNSS-PWV與ERA-PWV繪制散點圖,得到均方根誤差為2.28 mm,相關(guān)系數(shù)為0.910 3。這說明二者的擬合程度相對較好,證實了GNSS-PWV的準(zhǔn)確性和可靠性。

3)利用GNSS-PWV繪制大氣水汽含量變化圖,其結(jié)果的可信度表明,河南地區(qū)GNSS測站觀測數(shù)據(jù)可用于實時全天候的氣象預(yù)警工作中,以避免因出現(xiàn)此類特大暴雨造成的災(zāi)害和損失,同時也檢驗了GNSS測站在高時空分辨率特征下的良好性能。

4)暴雨發(fā)生期間一周內(nèi)的大氣可降水含量變化呈先增加后降低的趨勢,并且一直處于較高的狀態(tài)。結(jié)合氣象資料發(fā)現(xiàn),這主要受“煙花”和“查帕卡”臺風(fēng)攜大量海洋水汽、空氣不規(guī)則流動、河南地勢地形的影響。隨著GNSS氣象學(xué)的發(fā)展,GNSS-PWV可用于提高天氣預(yù)報的準(zhǔn)確度,也將有助于氣候變化研究和數(shù)值模型改進(jìn)。