聯(lián)合GNSS/RS多源數(shù)據(jù)反演三維大氣水汽分布研究

張文淵， 張書畢， 鄭南山， 張秋昭*， 丁楠

1 中國礦業(yè)大學(xué)自然資源部國土環(huán)境與災(zāi)害監(jiān)測重點(diǎn)實驗室， 江蘇徐州 221116 2 中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院， 江蘇徐州 221116 3 蘇黎世聯(lián)邦理工學(xué)院大地測量與攝影測量研究所， 蘇黎世 8093 4 江蘇師范大學(xué)地理測繪與城鄉(xiāng)規(guī)劃學(xué)院， 江蘇徐州 221116

0 引言

水汽是地球大氣的重要組成部分，是引起天氣變化的重要因素之一.水汽的相變會吸收或釋放巨大熱量，進(jìn)而引發(fā)臺風(fēng)、暴雨、干旱、洪水等破壞性的天氣現(xiàn)象.因此，準(zhǔn)確地獲取大氣水汽的時空信息對于改善天氣預(yù)報精度和災(zāi)害性氣象預(yù)警具有重要作用(Zhang et al., 2015; Wang et al., 2018; He et al., 2019).但由于水汽時空變化極為復(fù)雜，傳統(tǒng)探測手段受觀測成本、采樣間隔、測站分布等因素影響，難以實現(xiàn)對大氣水汽的實時高分辨率監(jiān)測(Zhang et al., 2015；姚宜斌等，2017b).Flores等(2000)首次提出GPS (Global Positioning System)水汽層析技術(shù)，利用夏威夷地區(qū)的GPS站網(wǎng)反演出局域三維濕折射率場.依靠其全天候監(jiān)測、高時空分辨率、低成本等優(yōu)勢，GNSS(Global Navigation Satellite System)層析技術(shù)在過去20年得到了快速發(fā)展(Champollion et al., 2005; 夏朋飛等, 2015; Benevides et al., 2017; 趙慶志等, 2018; Trzcina and Rohm, 2019; 張文淵等, 2021).三維濕折射率(wet refractivity, WR)和水汽密度(water vapor density, WVD)是層析技術(shù)的主要產(chǎn)品，目前已應(yīng)用于改善數(shù)值預(yù)報模式(Trzcina and Rohm, 2019)，研究臺風(fēng)周期變化(Zhu et al., 2020)，校正大氣延遲誤差(Yu et al., 2018)等領(lǐng)域.近年來，隨著GPS、BDS (BeiDou Navigation Satellite System)、GLONASS (GLObal NAvigation Satellite System)、Galileo (Galileo Satellite Navigation System)四大系統(tǒng)的發(fā)展和現(xiàn)代化，傳統(tǒng)GPS層析技術(shù)已發(fā)展為Multi-GNSS層析技術(shù)，利用多系統(tǒng)的GNSS觀測信號，可有效增加觀測值數(shù)量并改善GNSS信號的觀測幾何缺陷，提高層析結(jié)果的精度和可靠性(Benevides et al., 2017；Zhao et al., 2018；夏朋飛等，2015；胡鵬等，2020).

然而，受GNSS衛(wèi)星和測站空間分布的限制，GNSS信號穿刺的“倒錐形”觀測區(qū)域與“盒形”層析區(qū)域存在明顯的不匹配性，尤其是在近地層區(qū)域，導(dǎo)致大量的層析體素塊無法被GNSS信號線穿過，引起層析方程組的不適定性問題(Rohm, 2013)，這是GNSS水汽層析領(lǐng)域的研究難點(diǎn)和熱點(diǎn)(姚宜斌等，2017b).針對該問題，F(xiàn)lores等(2000)添加水平約束和垂直約束來估計空白體素塊的水汽參數(shù)；Song等(2006)利用高斯加權(quán)平均函數(shù)確定水平約束中周圍體素塊的權(quán)重，使得層析結(jié)果趨于真實的大氣水汽分布；何林等(2015)在水平方向采用均值濾波器進(jìn)行平滑處理，利用二階Laplacian算子確定周圍不同體素塊的權(quán)重信息；基于此，Zhang 等(2017)利用Helmert方差分量估計建立了自適應(yīng)Laplacian平滑方法確定最優(yōu)的權(quán)重，改進(jìn)后的層析結(jié)果與探空數(shù)據(jù)更為吻合.在層析模型優(yōu)化方面，Chen 和 Liu(2014)利用GNSS站網(wǎng)的分布特點(diǎn)，根據(jù)實驗確定最佳的層析模型水平分辨率，使得被穿刺的體素塊數(shù)量達(dá)到最大值；隨后，Rohm 等(2014)和Zhang 等(2015)共同提出內(nèi)外雙分辨率水汽層析方法，來降低層析區(qū)域邊界處空白體素塊的數(shù)量.Haji-Aghajany等(2020)利用WRF(Weather Research and Forecasting)模型獲得初始體素塊的先驗水汽參數(shù)，然后構(gòu)造新的優(yōu)化體素塊，大大地減少了體素塊的數(shù)量.姚宜斌等(2017a)和Brenot 等(2020)分別提出兩種附加虛擬觀測信號的精化水汽層析模型，通過增加穿過研究區(qū)域的射線條數(shù)和射線穿過的體素塊數(shù)，使層析結(jié)果更加逼近真值.

此外，融合外部氣象數(shù)據(jù)也是提高層析反演質(zhì)量的重要手段之一.Xia 等(2013)基于GPS觀測值和COSMIC(constellation observing system for meteorology, ionosphere and climate)廓線信息建立了兩步重構(gòu)法的GNSS水汽層析模型.Benevides等(2018)利用大氣紅外探測儀(Atmospheric Infrared Sounder, AIRS)獲取的高精度水汽廓線對GNSS層析模型進(jìn)行初始化，有效提高了層析解算結(jié)果的精度.趙慶志等(2018)和Zhang等(2020)分別利用ECMWF(European Centre for Medium-Range Weather Forecasts)的再分析資料和30年的Radiosonde數(shù)據(jù)構(gòu)造比例因子模型，用于精確估計側(cè)邊信號的水汽含量，大幅提高觀測信號數(shù)量和層析結(jié)果質(zhì)量.除此之外，高水平分辨率的遙感氣象數(shù)據(jù)也可以提供豐富的大氣水汽信息以及良好的幾何觀測結(jié)構(gòu)，例如，MODIS(Moderate-Resolution Imaging Spectroradiometer)、MERIS(Medium Resolution Imaging Spectrometer)、以及我國發(fā)射的FY-3A衛(wèi)星上的MERSI(Medium-Resolution Spectral Imager)傳感器都可以提供1 km或5 km水平分辨率的PWV觀測值(Chang et al., 2015; He and Liu, 2019)，這為融合RS和GNSS觀測數(shù)據(jù)提供了基礎(chǔ).

本文以聯(lián)合GNSS與RS多源觀測的水汽數(shù)據(jù)為切入點(diǎn)，充分利用兩種觀測信號的空間互補(bǔ)優(yōu)勢，構(gòu)造GNSS/RS聯(lián)合水汽層析觀測方程組，建立GNSS/RS聯(lián)合水汽層析模型，有效解決了傳統(tǒng)層析模型中觀測幾何限制的問題.利用徐州地區(qū)的GNSS實測數(shù)據(jù)和遙感水汽數(shù)據(jù)進(jìn)行實驗分析，并以研究區(qū)域內(nèi)的無線電探空數(shù)據(jù)和ERA5再分析資料對層析結(jié)果進(jìn)行精度評定，結(jié)果表明GNSS/RS層析模型可提高層析結(jié)果的精度，準(zhǔn)確獲取近地空間大氣水汽的三維分布情況.

1 水汽層析原理

水汽層析技術(shù)是指利用觀測信號在對流層中的濕延遲信息反演研究區(qū)域的三維濕折射率場或水汽密度場(Flores et al., 2000).在GNSS水汽層析中，觀測信號的斜路徑濕延遲(slant wet delay, SWD)可表示為濕折射率沿著信號方向從接收機(jī)到對流層頂?shù)姆e分：

(1)

式中，s表示GNSS信號的傳輸路徑，Nw為濕折射率.為反演研究區(qū)域的三維水汽密度分布，將SWD利用(2)式可轉(zhuǎn)換為斜路徑水汽含量(slant water vapor, SWV):

(2)

式中，105/RV·(k3/Tm+k′2)表示轉(zhuǎn)換因子，記為Π.其中k′2=16.48 K·hPa-1，k3=3.75×105K2·hPa-1，兩者均表示大氣折射率常數(shù)，RV表示水汽氣體常數(shù)，值為461.53 J/(kg·K) (Davis et al., 1985).Tm表示大氣加權(quán)平均溫度，本文利用Bevis模型計算Tm的值(Bevis et al., 1992)，該模型是目前應(yīng)用最為廣泛的Tm模型之一.根據(jù)水汽層析基本原理，將層析區(qū)域劃分為均勻分布的體素塊，則可建立每個體素塊的水汽參數(shù)與衛(wèi)星信號的SWV值的函數(shù)關(guān)系(Flores et al., 2000)：

(3)

式中，SWV單位為mm；dijk為斜路徑射線穿過位于第k層，第i行，第j列的體素塊的截距信息，單位為km；xijk表示層析模型的待求參數(shù)，即體素塊的水汽密度值，單位為g·m-3.由于層析模型中大量體素塊未被觀測信號穿過，這些體素塊所對應(yīng)的截距信息為0，導(dǎo)致了層析觀測方程的不適定性，通常采用水平和垂直經(jīng)驗函數(shù)對觀測方程進(jìn)行約束(Flores et al., 2000; Song et al., 2006; Benevides et al., 2017).

因此，聯(lián)合層析觀測方程、水平約束方程和垂直約束方程，可組成對流層層析方程組：

(4)

式中，SWVO和AO分別表示信號的SWV觀測向量和對應(yīng)的觀測矩陣；AH和AV表示水平約束和垂直約束的系數(shù)矩陣，X為水汽密度的未知參數(shù)向量.

2 GNSS/RS聯(lián)合水汽層析模型

2.1 GNSS/RS SWV聯(lián)合觀測值

由式(4)可知，精確的SWV觀測值是反演大氣水汽密度參數(shù)X的必要條件，GNSS/RS水汽層析模型的關(guān)鍵在于獲取高精度的GNSS SWV和RS SWV觀測值.圖1給出了GNSS/RS聯(lián)合水汽層析模型的具體流程，本節(jié)將詳細(xì)介紹這兩類關(guān)鍵觀測值的計算過程.

圖1 GNSS/RS聯(lián)合水汽層析模型流程圖Fig.1 Flow chart of the GNSS/RS combined water vapor tomography model

在GNSS觀測中，測站可以接收到來自不同衛(wèi)星的電磁波信號，利用雙差網(wǎng)解法或非差精密單點(diǎn)定位方法可以對GNSS測站處的對流層延遲和梯度信息進(jìn)行估計(范世杰, 2013)，然后利用(5)式可計算出觀測信號的SWD值(Flores et al., 2000)：

SWD=mfw(ele)·ZWD+mfg(ele)

(5)

在遙感觀測中，搭載在衛(wèi)星上的傳感器可以直接探測近地層的高分辨率水汽信息.例如，Terra和Aqua衛(wèi)星上的MODIS傳感器，Envisat衛(wèi)星上的MERIS傳感器，以及我國發(fā)射的FY-3A衛(wèi)星上的MERSI傳感器都可以提供高分辨率的PWV觀測值，并且具有全球觀測的優(yōu)勢(Chang et al., 2015; He et al., 2019).如圖2所示，與GNSS SWV的計算方法相似，可以利用下式估計RS像元到遙感衛(wèi)星的信號路徑上的RS SWV觀測值：

SWVRS(ε,α)=mw(ε)·PWVRS+Π·mfg(ε)

(6)

圖2 反演RS SWV觀測值的示意圖Fig.2 Schematic of the retrieval of RS SWV observations

由式(6)可知，獲取可靠的RS SWV觀測值的關(guān)鍵是確定精確的RS PWV和水平梯度延遲量.但是，由于傳感器的光譜信號易受云的影響，原始的RS PWV水汽產(chǎn)品存在觀測誤差.將其用于獲取RS SWV觀測值前，需要先對其進(jìn)行校正(Li et al., 2003; Chang et al., 2015).另一方面，大氣水汽的各向異性分布特點(diǎn)也會對斜路徑濕延遲量造成一定影響，本文選用GPT3w(Global Pressure and Temperature 3 wet)模型獲取大氣水平梯度分布信息(Landskron and B?hm, 2017).接下來將詳細(xì)論述RS PWV校正和梯度延遲估計的過程.

2.1.1 RS PWV校正

遙感水汽數(shù)據(jù)可以提供高分辨率的PWV觀測影像，以MODIS數(shù)據(jù)為例，其包括MOD05_L2和MYD05_L2兩種水汽產(chǎn)品，分別來自Terra和Aqua衛(wèi)星，均可以提供1 km水平分辨率的近紅外產(chǎn)品和5 km水平分辨率的紅外產(chǎn)品，通常近紅外數(shù)據(jù)只能在白天可以獲取，而紅外數(shù)據(jù)在白天和夜間均可以獲取，因此5 km的PWV數(shù)據(jù)的時間分辨率要高于1 km的PWV數(shù)據(jù).另外，GNSS層析模型的水平分辨率通常在10 km到50 km之間，若將1 km水平分辨率的PWV數(shù)據(jù)融合層析模型中，則會產(chǎn)生超大量的觀測信號.盡管增加了觀測值的數(shù)量，但是也會為層析方程解算帶來一定的負(fù)擔(dān)，而5 km的PWV數(shù)據(jù)既可以增加一定數(shù)量的觀測值，又不會加重計算負(fù)擔(dān).綜上，我們選用 5 km水平分辨率的PWV數(shù)據(jù).同時，為了保持一致性，需要將1 km水平分辨率的MERSI PWV數(shù)據(jù)稀釋成5 km水平分辨率的PWV產(chǎn)品.但由于遙感水汽產(chǎn)品觀測精度不高，通常以觀測區(qū)域內(nèi)同步GNSS數(shù)據(jù)獲取的GNSS PWV作為參考值，對RS PWV進(jìn)行校正.

利用GNSS ZWD和式(2)中的比例因子∏的乘積，可以獲得全天候監(jiān)測的GNSS PWV觀測值.大量實驗表明與探空站提供的高精度PWV觀測值相比，GNSS PWV的RMSE(root mean square error)為1～2 mm，可用于RS PWV的校正(Li et al., 2003; Chang et al., 2019).由于GNSS測站與RS像元未完全重合，一般利用反距離加權(quán)算法插值計算出GNSS測站處的RS PWV觀測值，并將其作為待檢驗值，然后利用同一時刻的GNSS PWV對其進(jìn)行校正(Chang et al., 2019; He et al., 2019).此外，考慮到大氣水汽的快速變化特性，本文采用逐影像校正的方法依次對所選RS PWV觀測影像構(gòu)造線性回歸校正模型，如下式所示(Li et al., 2003)：

(7)

2.1.2 大氣水平梯度延遲

(8)

式中，A0表示平均值，A1和B1表示年周期系數(shù)，A2和B2表示半年周期系數(shù)，這5個系數(shù)都已事先計算好，以格網(wǎng)形式保存在文本文件中(姚宜斌等, 2015)，doy表示年積日，需要結(jié)合RS影像的UTC時刻換算成精確的儒略日時間信息再進(jìn)行梯度估計.本文基于1°×1°的格網(wǎng)數(shù)據(jù)利用雙線性內(nèi)插得出RS影像獲取時刻每個RS像元的兩個梯度延遲量(朱明晨等, 2019).

2.2 GNSS/RS SWV聯(lián)合層析方程組

基于校正后的RS PWV觀測值和GPT3w模型估計的水平梯度延遲量，結(jié)合RS像元到衛(wèi)星的信號的高度角和方位角信息，利用式(6)可以準(zhǔn)確估計出RS信號的斜路徑水汽含量(RS SWV)，因此融合GNSS SWV和RS SWV的層析觀測方程如下所示：

(9)

式中，SWVGNSS和SWVRS分別表示GNSS SWV和RS SWV的觀測向量，AGNSS和ARS則表示由GNSS和RS信號在三維層析模型中的截距信息組成的觀測矩陣，AH和AV和與式(4)中一致.針對不同類型的觀測信號，通常需要對其進(jìn)行合理定權(quán)，參考Zhao 等(2017)和Heublein 等(2019)提出的定權(quán)方案，對GNSS和RS信號采用高度角定權(quán)，如式(10)所示；對于約束條件采用單位陣定權(quán).

Pi=sin2(ei),

(10)

式中，ei和Pi分別表示第i條觀測信號的高度角和權(quán)重信息.

代數(shù)重構(gòu)算法(Algebraic Reconstruction Technique, ART)采用逐行迭代的方式對線性方程組進(jìn)行解算，避免系數(shù)矩陣的求逆運(yùn)算，具有穩(wěn)定、高效的特點(diǎn)，在稀疏型方程解算中具有優(yōu)勢(Bender et al., 2011; 姚宜斌等, 2014; 張文淵等, 2021)，本文利用該算法對上述層析方程組進(jìn)行解算.

(11)

3 層析試驗分析

3.1 實驗數(shù)據(jù)

本文的實驗數(shù)據(jù)包括GNSS數(shù)據(jù)、MODIS和MERSI水汽數(shù)據(jù)、探空數(shù)據(jù)以及ERA5再分析數(shù)據(jù).研究時間為2017年8月，徐州地區(qū)正值夏季，多雷暴雨天氣.GNSS數(shù)據(jù)包括徐州Continuously Operating Reference Stations(CORS)網(wǎng)中10個均勻分布的GNSS測站提供的單系統(tǒng)GPS的觀測數(shù)據(jù)，GNSS測站的水平分布如圖3所示，其高程范圍為31.2 m到60.7 m，高程變化對三維層析模型的影響可忽略不計.利用GAMIT10.7軟件處理GNSS觀測數(shù)據(jù)，并引入3個IGS站(BJFS站，SHAO站，LHAZ站)的觀測數(shù)據(jù)，采用雙差解算策略進(jìn)行處理，以獲取高精度的ZTD觀測值(張文淵等, 2021).同時，采用Saastamoinen模型估計ZHD值(Saastamoinen, 1972)，選用VMF1濕映射函數(shù)計算GNSS SWV值(B?ehm and Schuh, 2004).利用GNSS數(shù)據(jù)估計ZTD和梯度信息的時間分辨率為5 min.

關(guān)于遙感水汽數(shù)據(jù)，本文采用Aqua、Terra兩個衛(wèi)星提供的MODIS數(shù)據(jù)和FY-3A衛(wèi)星提供的MERSI數(shù)據(jù).由于近幾年已不提供MERIS數(shù)據(jù)，所以文中沒有使用該數(shù)據(jù)(He et al., 2019).根據(jù)層析區(qū)域的觀測條件確定可用的MODIS和MERSI水汽數(shù)據(jù)，本文挑選置信度大于95%的RS PWV影像參與到層析實驗中(Li et al., 2003)，表1列出了符合要求的40幅影像的觀測衛(wèi)星和觀測時間.

圖3 徐州地區(qū)GNSS測站分布示意圖Fig.3 Distribution of GNSS stations in Xuzhou

表1 2017年8月徐州地區(qū)40幅RS PWV影像的觀測衛(wèi)星及獲取時間Table 1 Observation satellites and acquisition time of the 40 RS PWV maps in Xuzhou in August 2017

此外，本文采用探空數(shù)據(jù)和ERA5數(shù)據(jù)對層析結(jié)果進(jìn)行精度驗證.探空數(shù)據(jù)由美國俄懷明大學(xué)(University of Wyoming)提供，其包括不同高度層的氣壓、溫度、相對濕度、比濕等氣象數(shù)據(jù)，可以精確地計算出不同高度處的水汽密度含量，常用于驗證GNSS反演結(jié)果的精度(Li et al., 2003).ERA5數(shù)據(jù)是第五代大氣再分析資料，可以提供0.25°×0.25°水平分辨率的37個氣壓層逐小時分辨率三維氣象參數(shù).相比于傳統(tǒng)的ERA-Interim產(chǎn)品，ERA5數(shù)據(jù)具有更高的精度和時間分辨率(Cong et al., 2018).因此，我們分別采用高精度的探空垂直廓線和ERA5三維水汽分布作為驗證數(shù)據(jù)來計算GNSS/RS層析模型的反演精度.

3.2 層析策略

在本文的層析實驗中，層析區(qū)域為：經(jīng)度116.52—117.96°E，緯度33.90—34.62°N.參考現(xiàn)有研究中有關(guān)層析格網(wǎng)劃分的設(shè)定(Zhang et al., 2015, 2021; Benevides et al., 2018; Heublein et al., 2019; Haji-Aghajany et al., 2020)，考慮到徐州地區(qū)GNSS站網(wǎng)的稀疏分布特性，本文在經(jīng)度方向和緯度方向的水平分辨率均設(shè)置為0.18°，對應(yīng)的體素塊個數(shù)為8個和4個.垂直方向采用非均勻劃分的形式，共劃分15層，1到10層每層間隔為0.4 km，11到15層的間隔分別為1.0 km、1.0 km、1.5 km、1.5 km、2.0 km，層析頂層高度為11 km，根據(jù)徐州地區(qū)10年的探空歷史數(shù)據(jù)確定(張文淵等, 2021).單個層析歷元的時長為30 min，并且假設(shè)在一個層析歷元內(nèi)體素塊的水汽密度不發(fā)生變化(Benevides et al., 2018).本文層析時間段的設(shè)置主要根據(jù)遙感影像的觀測時間而定，使得RS PWV影像觀測時間恰好位于層析時段的中間時刻(例如，Image02的觀測時間為UTC 18∶30，則該影像對應(yīng)的層析歷元為UTC 18∶15—18∶45)，保證GNSS數(shù)據(jù)與RS數(shù)據(jù)的時間一致性.由于兩種遙感數(shù)據(jù)探測徐州地區(qū)的過境時間不一致，本文沒有將兩種遙感數(shù)據(jù)整合到一個層析歷元中，這也是遙感觀測的主要缺陷，其時間分辨率對于實時水汽反演有一定的限制.

本文設(shè)計如下兩種實驗方案來進(jìn)行對比分析，由于Terra、Aqua和FY-3A衛(wèi)星在徐州地區(qū)的過境時間并非每天的UTC 00∶00和UTC 12∶00，本文需要采用三次樣條插值對探空數(shù)據(jù)進(jìn)行時間域的插值(Chang et al., 2019)，以獲得層析時刻的探空水汽廓線信息.然后利用探空插值數(shù)據(jù)和ERA5再分析資料分別對兩種層析結(jié)果進(jìn)行精度評定，計算層析結(jié)果的RMSE、STD(standard deviation)和MAE(mean absolute error)等精度指標(biāo)，下邊給出精度指標(biāo)的計算公式：

方案一：只采用GNSS觀測信號進(jìn)行水汽反演，以“GNSS”表示；

方案二：采用GNSS與RS觀測信號聯(lián)合進(jìn)行水汽反演，以“GNSS/RS”表示.

(12)

(13)

(14)

3.3 RS PWV校正分析

根據(jù)2.1.1節(jié)介紹，選用逐影像的校正方法依次利用GNSS PWV觀測值對每幅遙感影像的PWV反演值進(jìn)行校正，表2列出了校正前后多項精度指標(biāo)的統(tǒng)計結(jié)果，可以看出校正后的RS PWV的精度有了明顯提高，其平均RMSE值從8.06 mm降低到1.05 mm，與GNSS PWV的相關(guān)系數(shù)由0.56提高到0.93，可以滿足水汽層析的要求(Heublein et al., 2019).

表2 校正前后RS PWV與GNSS PWV的相關(guān)系數(shù)以及精度統(tǒng)計量(mm)Table 2 Correlation coefficient and precision statistics of MODIS PWV and GNSS PWV before and after correction (mm)

3.4 層析模型觀測幾何分析

在傳統(tǒng)GNSS水汽層析模型中，受GNSS衛(wèi)星星座和測站分布的限制，位于層析模型低層區(qū)域的大量體素塊無法被GNSS信號穿刺.但是，不同于GNSS信號的倒錐形觀測形狀，RS信號的觀測幾何為正錐形，兩者恰好可以相互彌補(bǔ)，達(dá)到最優(yōu)觀測結(jié)構(gòu).圖4對比了2017年8月22日UTC 03∶10—03∶40層析歷元下，傳統(tǒng)GNSS層析模型和GNSS/RS聯(lián)合層析模型中觀測信號的三維空間分布.此外，圖5對比了該歷元下兩種層析模型中三維體素塊的穿刺情況.

圖4 2017年8月22日UTC 03∶10—03∶40傳統(tǒng)GNSS層析模型(a)和GNSS/RS聯(lián)合層析模型(b)信號分布圖，藍(lán)色線代表GNSS信號，紅色線代表RS信號Fig.4 Signal distribution of traditional GNSS tomography model (a) and GNSS/RS combing tomography model (b) at UTC 03∶10—03∶40 on August 22, 2017. Blue lines represent GNSS signals while red lines represent RS signals

圖5 傳統(tǒng)GNSS層析模型(a)和GNSS/RS聯(lián)合層析模型(b)的三維體素塊穿刺對比Fig.5 Comparison of 3D crossed voxels between the traditional GNSS tomography model (a) and GNSS/RS tomography model (b)

由圖5可知，在傳統(tǒng)GNSS層析模型中，低層存在大量的空白體素塊，而RS信號可以通過幾乎所有的體素塊.經(jīng)統(tǒng)計可知，GNSS層析模型中1到5層的穿刺體素塊數(shù)目為56個，占前5層總體素塊數(shù)量的35%，GNSS/RS模型的穿刺體素塊數(shù)目為160，所有體素塊均被穿刺.進(jìn)一步地，圖6比較了各個層析歷元內(nèi)兩種層析模型的觀測信號數(shù)量以及穿刺體素塊數(shù)量，其中平均觀測信號數(shù)量由GNSS模型的479條增加到GNSS/RS模型的663條；相應(yīng)地，平均穿刺體素塊數(shù)量由297個提升到462個.綜上可以看出加入RS信號可以充分穿刺層析模型的體素塊，有效改善傳統(tǒng)層析模型中觀測信號與層析區(qū)域不匹配的缺陷.

圖6 不同層析歷元下兩種模型的觀測信號數(shù)量(a)以及穿刺體素塊數(shù)量(b)對比Fig.6 Comparison of the number of observed signals (a) and the number of crossed voxels (b) for the two models during different tomographic epochs

3.5 層析結(jié)果精度分析

本節(jié)比較了GNSS層析模型和GNSS/RS層析模型反演的水汽分布情況，以探空水汽廓線和ERA5再分析資料為參考值，對兩種模型獲取的層析結(jié)果進(jìn)行了精度評定.如上所述，由于時間的不一致性，本文引入三次樣條插值算法對探空數(shù)據(jù)進(jìn)行插值，用于確定層析時刻的探空站參考數(shù)據(jù).盡管該算法已廣泛應(yīng)用于水汽產(chǎn)品的時間域插值處理，但不同的天氣條件對探空插值結(jié)果的可靠性有一定的影響.參考Heublein等 (2019)中關(guān)于探空插值數(shù)據(jù)可靠性的檢驗方法，當(dāng)滿足以下兩個條件時，探空插值廓線可以被視為可靠參考值：(1)在連續(xù)兩個探空觀測時刻上，GNSS PWV和探空PWV有很好的一致性；(2)探空站在層析時刻的插值PWV與同一時刻的GNSS PWV數(shù)據(jù)非常接近.在這兩個條件中，以2 mm作為PWV誤差的臨界值.圖7分別展示了2017年8月GNSS PWV與原始探空PWV的差值以及逐個層析時刻GNSS PWV與探空插值PWV的誤差，可以看出歷元1、3、4、12、14、15、16、17、22、25、26、27、28、29、35、36這16個層析時刻的探空插值數(shù)據(jù)符合要求.

圖7 2017年8月GNSS PWV與探空PWV的誤差(a)以及不同層析歷元GNSS PWV與探空插值PWV的誤差(b)Fig.7 (a) Difference between GNSS PWV and radiosonde PWV in August 2017 and (b) error of GNSS PWV and radiosonde interpolation PWV in different tomographic epochs

為了比較兩種水汽層析模型的層析結(jié)果的質(zhì)量，圖8首先展示了不同天氣條件下(晴天：2017年8月15日UTC 04∶40—05∶10，雨天：2017年8月18日UTC 22∶45—23∶15)兩種層析模型反演的大氣水汽廓線以及探空插值廓線(對應(yīng)上述兩個時間段).需要注意的是，由于探空站僅能提供其所在位置上空的垂直水汽廓線信息，因此這里只將探空站所處那一列的體素塊的層析水汽廓線與探空水汽廓線進(jìn)行對比.盡管通過驗證分析，探空插值數(shù)據(jù)仍會存在一定的誤差，不利用準(zhǔn)確分析層析結(jié)果的精度，因此引入獨(dú)立的ERA5數(shù)據(jù)進(jìn)行精度驗證，該數(shù)據(jù)還可以驗證層析結(jié)果的整體三維分布質(zhì)量，圖9 比較了相同時間段內(nèi)層析結(jié)果和ERA5水汽產(chǎn)品的三維分布情況.綜合圖8和9可以發(fā)現(xiàn)，在5 km(第11層)以上，由于兩種模型均有大量觀測信號穿過且分布較為均勻，對應(yīng)的兩種層析廓線較為接近，且與探空水汽廓線變化趨勢一致.但在5 km以下，GNSS層析模型反演的水汽廓線明顯低于GNSS/RS模型的水汽廓線，并且與探空水汽廓線也存在較大的差異.這可能是由于在GNSS層析模型中，“倒錐形”GNSS信號簇與“盒形”層析區(qū)域的空間幾何形態(tài)的不匹配性導(dǎo)致GNSS層析模型中大量的體素塊無法被GNSS信號線穿過.因此，在GNSS層析方程的解算過程中，這些空白體素塊的水汽信息只能通過約束方程來進(jìn)行估計，導(dǎo)致低層體素塊的水汽結(jié)果低于GNSS/RS模型和ERA5數(shù)據(jù).但在GNSS/RS層析模型中，通過引入“正錐形”RS信號觀測值，使得低層的絕大多數(shù)體素塊可以被信號穿過，保證三維層析模型的結(jié)構(gòu)穩(wěn)定性，進(jìn)而提高了層析結(jié)果的質(zhì)量.就層析模型的三維反演精度而言，同樣可以看出GNSS/RS層析模型得到的水汽三維分布與ERA5資料提供的水汽三維分布更加接近，進(jìn)一步表明遙感水汽信號對于提高三維水汽分布的反演精度有重要作用.

表3 兩種方法的層析結(jié)果的統(tǒng)計量對比 (單位：g·m-3)Table 3 Statistics of the tomography results for two methods (unit: g·m-3)

圖8 兩種層析模型在不同天氣條件下的水汽廓線對比(a) Rainless天氣； (b) Rainy天氣.Fig.8 Comparison of water vapor profiles for two tomography models under different weather conditions(a) Rainless weather; (b) Rainy weather.

圖9 兩種層析模型在不同天氣條件下的三維水汽分布對比Fig.9 Comparison of 3D water vapor distribution for two tomography models under different weather conditions

為了更好的比較不同層析歷元下兩種模型的反演精度，圖10給出了層析時段內(nèi)兩種模型的層析結(jié)果與探空數(shù)據(jù)和ERA5數(shù)據(jù)的RMSE對比.由圖可知，在絕大多數(shù)層析歷元內(nèi)，GNSS/RS層析模型的RMSE都明顯優(yōu)于GNSS模型.此外，表3列出了不同模型的層析結(jié)果的精度統(tǒng)計量，可以看出以探空站為參考值，GNSS模型和GNSS/RS模型的平均RMSE分別為2.19 g·m-3、1.54 g·m-3，平均MAE分別為1.69 g·m-3、1.34 g·m-3，兩者分別提高了29.68%和20.71%；以ERA5數(shù)據(jù)為參考值評估40個歷元的層析結(jié)果可得兩個模型的平均RMSE分別為2.56 g·m-3、2.04 g·m-3，平均MAE分別為2.03 g·m-3、1.63 g·m-3，分別提高了20.31%和21.18%.綜合來看，相對于傳統(tǒng)GNSS層析模型，本文提出的GNSS/RS層析模型的RMSE和MAE分別提高了25.00%和20.95%，這表明融合RS信號可以有效提高層析結(jié)果的精度.

圖10 層析時段內(nèi)兩種模型的層析結(jié)果分別與有效探空插值數(shù)據(jù)(a)和ERA5數(shù)據(jù)(b)對比的RMSE值Fig.10 RMSE values of the tomographic results from the two models are compared with effective radiosonde interpolation data (a) and ERA5 data (b) respectively in the tomographic periods

為了進(jìn)一步比較不同層析方法在垂直方向上的反演精度，分別以探空插值數(shù)據(jù)和ERA5數(shù)據(jù)為參考值，對不同高度層上的層析結(jié)果進(jìn)行對比分析.圖11給出了不同高度層上，兩種層析水汽結(jié)果與探空數(shù)據(jù)和ERA5數(shù)據(jù)對比的均方根誤差(RMSE)、平均偏差(mean deviation, MD)、標(biāo)準(zhǔn)偏差(STD)以及相對誤差(relative error, RE)四個統(tǒng)計指標(biāo)信息.可以發(fā)現(xiàn)，相對于兩種參考數(shù)據(jù)，GNSS/RS聯(lián)合模型的層析結(jié)果質(zhì)量在絕大多數(shù)高度層上都優(yōu)于傳統(tǒng)GNSS模型，尤其是在0～2 km的近地范圍內(nèi)，改善效果最為明顯.該范圍內(nèi)的層析結(jié)果的平均RMSE值由2.95 g·m-3降低到2.04 g·m-3，平均偏差和標(biāo)準(zhǔn)偏差分別減小了1.11 g·m-3和0.75 g·m-3，相對誤差由17.45%改善為10.84%.這進(jìn)一步證明了GNSS/RS層析方法在反演高精度三維水汽分布方面的優(yōu)勢.

圖11 層析時段內(nèi)兩種模型的層析結(jié)果分別與探空插值數(shù)據(jù)(a—d)和ERA5數(shù)據(jù)(e—h)對比的分層精度信息Fig.11 Layered statistics of the tomographic results from the two models compared with radiosonde interpolation data (a—d) and ERA5 data (e—h) in the tomographic periods

4 結(jié)論

本文基于高水平分辨率RS信號的觀測優(yōu)勢，首次創(chuàng)建了聯(lián)合GNSS與RS多源觀測的水汽層析模型，并詳細(xì)論述了GNSS/RS層析模型的構(gòu)造流程.利用2017年8月徐州地區(qū)的GNSS實測數(shù)據(jù)和MODIS、MERSI水汽遙感數(shù)據(jù)進(jìn)行水汽層析實驗，系統(tǒng)地比較分析了傳統(tǒng)GNSS模型與GNSS/RS模型的觀測幾何結(jié)構(gòu)和層析結(jié)果質(zhì)量，主要結(jié)論如下：

(1)RS PWV作為獲得準(zhǔn)確RS SWV觀測值的前提條件，在進(jìn)行層析之前必須對其進(jìn)行校正.利用同步的GNSS PWV對其進(jìn)行校正后，RS PWV的平均RMSE由8.06 mm降低為1.05 mm，可以滿足層析試驗的要求.

(2)在層析模型的觀測幾何結(jié)構(gòu)方面，GNSS/RS層析模型的平均觀測信號數(shù)量和穿刺體素塊數(shù)量分別為663條和462個，明顯優(yōu)于傳統(tǒng)GNSS模型的479條和297個，兩者分別提高了38.41%和55.56%.在GNSS/RS模型中，位于低層的多數(shù)體素塊均能被RS信號較好地穿刺，有效改善傳統(tǒng)層析模型的GNSS信號幾何結(jié)構(gòu)缺陷.

(3)分別以探空水汽廓線和ERA5三維水汽分布信息為參考值對兩種模型的層析結(jié)果進(jìn)行評估.

發(fā)現(xiàn)層析水汽廓線的平均RMSE從GNSS模型的2.19 g·m-3降低到GNSS/RS模型的1.54 g·m-3，精度提高了29.68%；就層析三維水汽分布情況而言，傳統(tǒng)模型的RMSE為2.56 g·m-3，多源觀測模型為2.04 g·m-3，三維分布精度改善了20.31%.此外，在不同高度層上，GNSS/RS聯(lián)合模型的層析結(jié)果的精度統(tǒng)計量也優(yōu)于傳統(tǒng)的GNSS層析模型.結(jié)合本文可知高分辨率遙感數(shù)據(jù)在改善層析反演結(jié)果方面展現(xiàn)出較大的潛力，受其時間分辨率的影響，目前對于實時水汽反演存在一定的局限性.但是，遙感數(shù)據(jù)的大范圍、高空間分辨率的優(yōu)勢，有助于進(jìn)行大尺度水汽層析實驗和三維水汽場精細(xì)化建模，這將是我們下一步的重點(diǎn)研究工作.

致謝作者特別感謝徐州市自然資源和規(guī)劃局提供的XZCORS網(wǎng)監(jiān)測數(shù)據(jù)，美國國家航空航天局提供的MODIS PWV影像，中國風(fēng)云衛(wèi)星遙感數(shù)據(jù)服務(wù)網(wǎng)提供的MERSI PWV影像，美國俄懷明大學(xué)提供的長期探空數(shù)據(jù)，麻省理工學(xué)院提供的GAMIT/GLOBK軟件.