聯(lián)合地基GNSS及空基GNSS掩星探測水汽三維分布

吳 昊，鄂盛龍，夏朋飛，許海林，羅穎婷，譚理慶

(1.廣東電力科學研究院，廣州 510080；2.武漢大學GNSS技術(shù)研究中心，武漢 430072)

0 引言

全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satel-lite System，GNSS)能夠有效地服務(wù)于全球?qū)Ш脚c定位，此外，GNSS技術(shù)又擴展應(yīng)用到全球氣候監(jiān)測與極端天氣的預(yù)報[ 1-3]。目前，利用GNSS技術(shù)探測地球大氣已成為國內(nèi)研究的熱點之一[ 1-5]。作為一種新的大氣探測手段，GNSS技術(shù)能夠準實時、無需人為干擾、全天候地服務(wù)大氣科學的研究[ 6-10]。根據(jù)觀測模式的不同，GNSS探測技術(shù)可以分為空基GNSS技術(shù)和地基GNSS技術(shù)[ 8-10]?？栈鵊NSS掩星技術(shù)能夠提供高垂直分辨率、高精度、全球覆蓋的中性大氣層到電離層的多級大氣剖面產(chǎn)品。然而，受制于衛(wèi)星星座的分布及數(shù)量，目前空基掩星技術(shù)還無法實時地服務(wù)于全球數(shù)值天氣預(yù)報模式[10]。

地基GNSS系統(tǒng)探測大氣是利用平差的方法從觀測數(shù)據(jù)中估算出天頂對流層延遲(Zenith Tropospheric Delay，ZTD)，將其作為研究對流層大氣的背景場。通常ZTD可以分為2個部分，即天頂對流層干延遲(Zenith Hydrostatic Delay，ZHD)和天頂對流層濕度延遲(Zenith Wet Delay，ZWD)。ZWD中含有大氣水汽的信息，因此，利用一定的數(shù)學模型可以從ZWD中反演出大氣水汽的分布[ 8-10]。

為了提高GNSS技術(shù)探測大氣水汽分布的質(zhì)量，文中試圖利用無線電掩星產(chǎn)品輔助地基GNSS技術(shù)反演水汽三維分布信息。對流層天頂干延遲模型、濕度轉(zhuǎn)換因子以及層析模型上邊界的高度是影響層析結(jié)果的3個關(guān)鍵的參數(shù)。首先，根據(jù)無線電掩星提供的高垂直分辨率的大氣產(chǎn)品標定了對流層干延遲模型；然后，建立新的區(qū)域性大氣加權(quán)平均溫度模型，提高了濕度轉(zhuǎn)換因子的精度；最后，根據(jù)掩星提供的大氣產(chǎn)品確定了水汽層層頂，將其作為層析模型的上邊界層。利用香港地區(qū)2017年6月份的數(shù)據(jù)進行試算和分析，反演結(jié)果驗證了上述方法的可行性。

1 標定對流層干項延遲及濕度轉(zhuǎn)換因子

與對流層濕項延遲相比，對流層干項變化較為平穩(wěn)，通過建立經(jīng)驗?zāi)Ｐ涂梢暂^為精確的估算。國內(nèi)外研究者建立了許多全球或區(qū)域性的天頂對流層干延遲經(jīng)驗?zāi)Ｐ?。其中，較為常用的主要包括如下3個：Saastamoinen模型、Hopfield模型和Black模型[ 11-13]。然而，經(jīng)典的對流層干延遲模型存在著系統(tǒng)性的偏差，可以利用高精度的氣象產(chǎn)品對其進行修正，如探空和掩星產(chǎn)品。利用高精度氣象產(chǎn)品估算ZHD的數(shù)學公式如下[9]

(1)

(2)

其中，k1為常數(shù)，其經(jīng)驗值通常取為k1=77.6K/mbar；Pd為大氣干壓強；Pw為大氣水汽壓；T為溫度；Rd和Rw分別為干普適氣體常數(shù)和濕普適氣體常數(shù)。

選取2007—2016年中六月份與香港地區(qū)共址的COSMIC (The Constellation Observing System for Meteorology, Ionosphere, and Climate)和GRACE (Gravity Recovery and Climate Experiment)掩星事件，將wetPrf產(chǎn)品利用式(1)和式(2)的計算結(jié)果分別與三種經(jīng)典ZHD模型的計算結(jié)果進行了比較，其統(tǒng)計結(jié)果如表1所示。

表1 掩星產(chǎn)品計算的ZHD與經(jīng)典ZHD模型計算的ZHD的差Tab.1 Bias of ZHD between radio occultation-derived and ZHD models-derived cm

從表1可以看出，將無線電掩星產(chǎn)品估算的ZHD(ROZHD) 視為參考值，三種經(jīng)驗?zāi)Ｐ团c參考值之間都存在一定的偏差。其中，Sasstamoinen模型估算的ZHD (SAAS) 與參考值之間的偏差優(yōu)于0.93cm，這個精度遠遠高于Hopfield模型和Black模型估算的ZHD(HOPE和BLA)。此外，HOPF的精度要優(yōu)于BLA。為了提高ZHD的估算精度，文中試圖利用Hopfield模型的表達形式來修正Sasstamoinen模型。

(3)

這里，a、b和c為區(qū)域性標定模型的系數(shù)。通過2007—2016年6月份與香港地區(qū)共址的COSMIC和GRACE掩星產(chǎn)品來估算這3個參數(shù)。精確獲取ZHD后，便可從ZTD中分離出ZWD，利用濕度轉(zhuǎn)換因子便可得到大氣水汽含量，而濕度轉(zhuǎn)換因子是關(guān)于大氣加權(quán)平均溫度的函數(shù)[14-15]

(4)

(5)

其中，ρv為水的密度，且ρv=103kg/m3；k1、k2、k3都為常數(shù)，取值分別為k1=77.6K/hPa，k2=70.4K/hPa，k3=3.739×105K/hPa；md和mw分別為干和濕大氣的摩爾質(zhì)量，取值分別為md=28.96kg/kmol，mw=18.02kg/kmol；R表示普適氣體常數(shù)，取值為R= 8314Pa·m3·K-1·kmol-1。

根據(jù)所選取的香港地區(qū)2007—2016年的COSMIC和GRACE無線電掩星產(chǎn)品(RO)，利用式(5)便可估算出大氣加權(quán)平均溫度Tm，建立二階擬合關(guān)系如圖1所示。

圖1 大氣加權(quán)平均溫度二階擬合圖Fig.1 Second-order fitting of atmospheric weighted average temperature

從圖1可以看出，利用二階模型擬合的大氣加權(quán)平均溫度模型與掩星產(chǎn)品獲得的Tm具有較好的一致性，二者偏差的RMS值優(yōu)于2.43K。圖1擬合出的結(jié)果便為相關(guān)地區(qū)大氣加權(quán)平均溫度的模型

(6)

式中，Ts表示地表處的大氣溫度。

為了檢驗新的Tm模型的精度，文中將45004探空測站2017年的探空產(chǎn)品利用式(5)計算獲得的Tm作為真值。此外，文中又分別評估了Bevis模型[16](Tm=70.2+0.72Ts)及Wang[17]等香港區(qū)域性模型(Tm=113.29+0.5863Ts)計算的Tm與探空資料估算Tm的差值，如表2所示。

表2 模型估算的Tm與真值的差值Tab.2 Bias of Tm between models-derived and radiosonde-derived K

表2顯示了新的模型(NTm)估算的結(jié)果要明顯優(yōu)于Bevis模型結(jié)果(BTm)和Wang模型結(jié)果(WTm)。新的模型與真值之間的偏差的RMS值優(yōu)于2.94K，偏差的變化幅度為[-3.70, 7.19]。

2 新的確定水汽層層頂?shù)姆椒?/h2>

層析網(wǎng)格劃分的上下限的高度是從地面到對流層頂。實際上，從某個高度到對流層頂?shù)膮^(qū)域中水汽密度接近于0[9]。由于其值太小，在實際的反演過程中得到的反演結(jié)果通常出現(xiàn)負值或異常值。為了提高對流層層析的精度，有學者定義了某個低于對流層頂?shù)母叨茸鳛閷游龅纳舷薷撸⒚麨樗麑訉禹擺9-10]，將其作為水汽反演的上限高度。為了精確地確定水汽層層頂?shù)母叨龋疚奶岢隽艘环N新的方法確定水汽層層頂。

在某個高度h上的大氣可降水量PWVh可以表示為水汽密度的積分

(7)

其中，ρw表示液態(tài)水密度；ρk表示水汽密度。

則在垂直方向上，PWV呈現(xiàn)近似高斯指數(shù)函數(shù)分布

(8)

其中，PWV0為地表處的PWV；H為大氣層等效高度。

通過對式(8)求導，便可得到不同高度處的斜率Kh

Kh=tan(α)

(9)

理論上，當Kh接近為0時，便對應(yīng)為水汽層層頂?shù)母叨?。為了確定水汽層層頂?shù)母叨?，需要找到合適的Kh。本文定義當α=179.9°時對應(yīng)的高度為水汽層層頂。通過香港地區(qū)2007—2016年6月份的掩星反演產(chǎn)品計算出香港地區(qū)水汽層層頂高度的變化，如圖2所示。

圖2 利用掩星產(chǎn)品定義的香港地區(qū)水汽層層頂高度隨時間變化的規(guī)律Fig.2 Determining the water vapor top using radio occultation products in Hong Kong

從圖2可以看出，香港地區(qū)6月份水汽層層頂變化范圍為9.0km～10.6km。為了方便計算，文中求取水汽層層頂?shù)钠骄?.5km作為香港地區(qū)6月份的水汽層層頂。

3 試驗分析

選取香港地區(qū)12個連續(xù)氣象監(jiān)測站2017年6月份的觀測數(shù)據(jù)進行試算。利用Bernese 5.2軟件進行處理獲取天頂對流層延遲ZTD，基于式(3)對Sasstamoinen模型進行改正后得到ZHD，然后從ZTD中分離出ZWD。根據(jù)式(6)估算的大氣加權(quán)平均溫度代入式(4)中來獲取濕度轉(zhuǎn)換因子。選用Niell投影函數(shù)將ZWD投影到斜路徑方向上SWD后，利用濕度轉(zhuǎn)化因子將SWD轉(zhuǎn)換為斜路徑水汽含量SWV。把研究的香港地區(qū)在水平方向上劃分5×5=25個網(wǎng)格，如圖3所示。在垂直方向上，0.8km以下每400m劃分一層，共2層；在0.8km～8.6km以上每600m劃分一層，共13層；8.6km～9.5km劃分一層，故垂直網(wǎng)格共劃分成16個網(wǎng)格，則總網(wǎng)格為16×5×5=400個。選取Flores提出的經(jīng)典的對流層層析方法，便可反演出水汽的三維分布信息。

圖3 水平方向?qū)游鼍W(wǎng)格劃分及測站的分布Fig.3 Horizontal distribution of GNSS stations and the division of voxels

根據(jù)上述的數(shù)據(jù)處理方式，每0.5h輸出一次結(jié)果。限于篇幅，圖4僅僅給出了第160d的處理結(jié)果。

圖4 水汽的三維分布信息Fig.4 3D distribution of water vapor

由圖4可以看出，2017年第160d的水汽主要集中在5km以下，而在5km以上水汽含量較少。隨著時間的變化，0.8km～5km處的水汽密度發(fā)生較為明顯的變化。為了檢驗優(yōu)化后的層析結(jié)果，文中又利用傳統(tǒng)的層析方法反演出對流層水汽密度值。將45004探測站的觀測產(chǎn)品作為真值分別檢驗優(yōu)化后的層析結(jié)果(Opti)和傳統(tǒng)的層析結(jié)果(Trad)，如圖5所示。

圖5 層析結(jié)果與探空產(chǎn)品比較Fig.5 Statistical results between Opti & Trad and Radisonde

由圖5可以看出，層析反演的水汽密度值與探空產(chǎn)品的變化趨于一致。當水汽密度在空間上存在逆增層的情況下，層析水汽技術(shù)不能很好地反演出這種情況；Opti結(jié)果要稍好于Trad。此外，從圖5中還可以看出，Opti、Trad與Radiosonde在5km以上符合性較好，在5km以下差距較大。為了統(tǒng)計層析技術(shù)估算的水汽密度與探空水汽密度值之間的差距，文中分段分別統(tǒng)計了2017年6月份5km以下和5km以上的層析水汽值與探空水汽值的差值，如表3所示。

表3 Opti和Trad與Radiosonde偏差統(tǒng)計圖Tab.3 The statistical results between Opti & Trad and Radisonde

由表5可以看出，在5km以下由于存在逆增層現(xiàn)象，故Opti、Trad與Radiosonde的偏差較大，而在5km以上偏差較小。Opti與Trad相比，在5km以下精度提高了1.24g/m3，在5km以上提高了0.22g/m3。

4 結(jié)論

利用無線電掩星產(chǎn)品輔助地基GNSS探測大氣水汽的空間分布。針對對流層層析過程中3個關(guān)鍵技術(shù)進行優(yōu)化：1)利用無線電掩星歷史資料標定傳統(tǒng)的Sasstamoinen干延遲模型，并建立了區(qū)域性大氣加權(quán)平均溫度模型，提高濕度轉(zhuǎn)換因子的精度；2)精確地估算了水汽層層頂?shù)母叨?，將其作為層析網(wǎng)格的上邊界層。選取了香港地區(qū)12個監(jiān)測站2017年6月份的觀測數(shù)據(jù)進行試算，將反演的結(jié)果與探空產(chǎn)品進行比較和統(tǒng)計，驗證了優(yōu)化后層析方法的可行性和有效性。

致謝非常感謝香港特別行政區(qū)地政署測繪局提供的GNSS觀測數(shù)據(jù)及氣象觀測數(shù)據(jù)。非常感謝UCAR提供的COSMIC和GRACE掩星觀測產(chǎn)品。