GNSS/MODIS信號緊耦合水汽層析算法

張文淵，張書畢，鄭南山，丁 楠，劉 鑫，馬朋序

1. 中國礦業(yè)大學(xué)自然資源部國土環(huán)境與災(zāi)害監(jiān)測重點實驗室，江蘇 徐州 221116； 2. 中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇 徐州 221116； 3. 江蘇師范大學(xué)地理測繪與城鄉(xiāng)規(guī)劃學(xué)院，江蘇 徐州 221116

水汽是大氣的重要組成成分，具有極為復(fù)雜的時空變化特性，在全球氣候變化以及極端氣象災(zāi)害中都扮演著重要角色。地基GNSS水汽層析技術(shù)憑借其全天候監(jiān)測、高時空分辨率、不受天氣影響以及監(jiān)測成本低等優(yōu)勢，已成為探測大氣水汽高時空分辨率變化特性的重要手段之一[1]。

GNSS水汽層析技術(shù)是指利用區(qū)域地基CORS(continuously operating reference stations)網(wǎng)進行局域?qū)α鲗尤S水汽場的反演，由文獻[2]首次提出，并重構(gòu)出夏威夷地區(qū)的三維濕折射率分布。在傳統(tǒng)GNSS水汽層析模型中，由于“盒形”層析區(qū)域與“倒錐形”GNSS信號簇的空間幾何形態(tài)的不匹配性，導(dǎo)致三維層析模型中大量的體素塊無法被GNSS信號線穿刺，進而引起層析觀測方程的秩虧性，這是GNSS水汽層析領(lǐng)域的研究難點和熱點[1]。針對該問題，眾多學(xué)者提出了不同的優(yōu)化模型與算法[2-19]。文獻[2]充分考慮水汽的空間分布特性，構(gòu)造了水平約束和垂直約束條件來對空白體素塊進行約束處理。文獻[3—4]分別基于高斯距離加權(quán)函數(shù)與圖像平滑原理對水平約束條件優(yōu)化，使得三維水汽分布更符合真實狀態(tài)。文獻[5—9]基于GPS、BDS、GLONASS、Galileo多系統(tǒng)的衛(wèi)星觀測信號用于三維水汽層析研究，通過增加GNSS觀測信號，大大降低了空白體素塊的數(shù)量，并提高層析結(jié)果的精度。文獻[10—11]基于虛擬思想，構(gòu)造出虛擬GNSS觀測站和觀測信號來增加觀測信號數(shù)量，以此來提高層析模型體素塊的穿刺率。文獻[12—14]對GNSS層析模型的空間分辨率進行深入研究，通過優(yōu)化水平分辨率與垂直分辨率使GNSS信號穿刺盡可能多的體素塊，改善層析結(jié)果的質(zhì)量。此外，融合外部氣象數(shù)據(jù)也是改善層析秩虧問題的重要手段之一。文獻[15]利用探空數(shù)據(jù)和大氣紅外探測儀(atmospheric infrared sounder,AIRS)獲取的高精度水汽廓線作為先驗信息來估計三維水汽分布。文獻[16]根據(jù)COSMIC(constellation observing system for meteorology,ionosphere and climate)提供的水汽信息構(gòu)造了基于兩步重構(gòu)法的GNSS層析模型，結(jié)合先驗水汽觀測值有效提高了層析解算結(jié)果的精度。文獻[17—18]利用ECMWF(European Centre for Medium-Range Weather Forecasts)的再分析資料和多年的Radiosonde數(shù)據(jù)構(gòu)造比例因子模型，用于精確估計側(cè)邊信號的水汽含量，大幅提高觀測信號數(shù)量和層析結(jié)果質(zhì)量。

近年來，隨著各類遙感衛(wèi)星的升空，大量的水汽遙感數(shù)據(jù)憑借其監(jiān)測范圍廣，空間分辨率高的優(yōu)勢已被廣泛用于大范圍對流層水汽變化監(jiān)測研究。文獻[20]利用大范圍的InSAR數(shù)據(jù)獲取高分辨率的可降水量(precipitable water vapor,PWV)差分觀測值，并通過試驗證明InSAR數(shù)據(jù)用于大氣監(jiān)測的有效性與可靠性，這為高分辨率遙感數(shù)據(jù)用于GNSS水汽層析模型研究提供了可能。文獻[21—22]利用上述的PWV差分觀測值構(gòu)造新型GNSS層析模型約束條件，提高GNSS層析結(jié)果的精度。但由于InSAR數(shù)據(jù)提供的是差分PWV信息，在與GNSS信號融合過程中需要將差分PWV轉(zhuǎn)換為完整的PWV數(shù)據(jù)，一般借助于數(shù)值天氣預(yù)報模型(numerical weather prediction,NWP)提供的全球水汽背景場進行轉(zhuǎn)換[23]，這會降低高分辨率PWV觀測值的精度。與InSAR數(shù)據(jù)類似的有MODIS、MERIS、FY-3A等水汽遙感數(shù)據(jù)，該類遙感數(shù)據(jù)可直接提供高分辨率的完整PWV觀測值[23-24]，可被視為融合高分辨率遙感數(shù)據(jù)的GNSS水汽層析技術(shù)研究的理想數(shù)據(jù)源。

本文基于高分辨率MODIS PWV觀測信息，充分挖掘其空間幾何結(jié)構(gòu)優(yōu)勢，構(gòu)造出GNSS/MODIS緊耦合水汽層析算法。該算法首先利用GNSS PWV對MODIS PWV進行檢驗和校正，并分析了傳統(tǒng)GNSS層析模型利用MODIS信號的缺點，提出基于體素節(jié)點參數(shù)化構(gòu)造GNSS/MODIS緊耦合層析觀測方程組的方法，可充分利用MODIS觀測信號，并有效緩解GNSS層析模型的病態(tài)性問題。此外，利用徐州地區(qū)的探空數(shù)據(jù)對層析結(jié)果進行了精度評估，結(jié)果表明本文提出的緊耦合算法可以增強GNSS層析模型的穩(wěn)定性，保證重構(gòu)三維水汽場的準(zhǔn)確性。

1 GNSS水汽層析原理

GNSS信號的斜路徑水汽總量(slant water vapor,SWV)分為各向同性與各向異性兩部分。利用濕映射函數(shù)VFM1[25]和水平濕延遲梯度可計算出信號的斜路徑濕延遲(slant wet delay,SWD)，再通過式(2)將其轉(zhuǎn)換為SWV

(1)

SWV=Π·SWD

(2)

(3)

在GNSS層析模型中，通過對三維層析區(qū)域進行離散化可以獲得均勻分布的體素塊，一般用體素塊中心位置處的水汽密度值代表該體素塊的水汽密度，并且假設(shè)在一個層析歷元內(nèi)體素塊的水汽密度為常數(shù)。因此衛(wèi)星信號路徑上的SWV可由式(4)表達

(4)

式中，SWV表示衛(wèi)星信號路徑上的水汽含量值，層析區(qū)域一共分成t層，每一層內(nèi)按一定的經(jīng)度差和緯度差均勻劃分為n×m個體素塊；dijk表示衛(wèi)星信號穿過位于第k層，第i行，第j列的體素塊的截距；xijk表示對應(yīng)體素塊的水汽密度。

2 MODIS水汽產(chǎn)品檢驗及校正

中分辨率成像光譜儀(moderate resolution imaging spectroradiometer,MODIS)是美國發(fā)射的兩顆環(huán)境遙感衛(wèi)星Terra和Aqua上搭載的主要探測儀器之一，可提供36個光譜波段的地球觀測信息，為自然災(zāi)害與生態(tài)環(huán)境監(jiān)測、全球環(huán)境和氣候變化研究提供強有力的數(shù)據(jù)支持。這兩顆衛(wèi)星兩天可對全球進行一次觀測，可滿足突發(fā)性、快速變化的環(huán)境監(jiān)測[28]。此外，憑借MODIS數(shù)據(jù)高時空分辨率的優(yōu)勢，其提供的高分辨率PWV數(shù)據(jù)已被廣泛應(yīng)用于InSAR數(shù)據(jù)的大氣延遲改正[28]。但由于MODIS水汽產(chǎn)品具有一定的滯后期，且時間分辨率明顯低于GNSS的全天候監(jiān)測，其在實時大氣水汽反演中存在一定的局限性。

2.1 MODIS水汽產(chǎn)品簡介

MODIS數(shù)據(jù)包括Terra衛(wèi)星平臺提供的MOD05_L2以及Aqua衛(wèi)星提供的MYD05_L2兩種水汽產(chǎn)品，包含每個像元的經(jīng)緯度、掃描時間、太陽天頂角和方位角、傳感器天頂角和方位角、近紅外水汽值、紅外水汽值、云掩膜產(chǎn)品等13個參數(shù)[24]。本文采用來自Terra平臺的MOD05_L2水汽產(chǎn)品進行層析試驗。

MOD05_L2產(chǎn)品數(shù)據(jù)的空間分辨率有兩種：1 km×1 km和5 km×5 km。在InSAR大氣延遲改正中，通常選用與SAR影像接近的1 km×1 km高分辨率水汽數(shù)據(jù)。在本文中，一方面考慮到層析模型的水平分辨率一般大于10 km×10 km[2-19]；另一方面，1 km×1 km分辨率的MODIS數(shù)據(jù)會成倍地增加GNSS層析模型的計算量，因此選用5 km×5 km分辨率的水汽產(chǎn)品用于層析試驗研究。圖1(a)和(b)分別顯示了2016年7月1日14:00時我國東部地區(qū)的MODIS水汽數(shù)據(jù)以及徐州地區(qū)的局部放大圖，圖1(a)中的灰色區(qū)域表示被云朵污染的部分。

2.2 MODIS PWV檢驗與校正

如圖1(a)所示，MODIS數(shù)據(jù)可以提供每個像素點的可降水汽含量，簡記為MODIS PWV。由于MODIS信號易受云朵的影響，需要利用MODIS的云掩膜產(chǎn)品將受污染的MODIS PWV剔除，此外，還需利用高精度GNSS PWV對有效的MODIS PWV觀測值進行檢驗和校正[29]。本文選用徐州地區(qū)的5個GNSS測站數(shù)據(jù)和對應(yīng)時間上通過該地區(qū)的MODIS PWV數(shù)據(jù)進行比較分析。GNSS測站分布如圖1(b)所示，ZHJI、LGUO、BNTG、SGOU、SANP代表5個測站的位置，RS 58027對應(yīng)的是探空站的位置分布。同時，選取徐州地區(qū)2016年7月中15幅觀測質(zhì)量較好的MODIS影像進行云朵影響剔除，然后參與層析試驗。表1列出了15幅MODIS影像的獲取時間和對應(yīng)的層析時段范圍。

圖1 我國東部地區(qū)及徐州地區(qū)MODIS水汽分布圖Fig.1 MODIS water vapor map of eastern China and Xuzhou, Jiangsu province

表1 2016年7月徐州地區(qū)15幅MODIS影像獲取時間及層析時段

基于GNSS全天候運行的優(yōu)勢，可以通過處理GNSS數(shù)據(jù)獲取與MODIS影像同一時刻的水汽信息，利用轉(zhuǎn)換系數(shù)可以將式(1)中的ZWD轉(zhuǎn)換為PWV，如式(5)所示

PWV=Π·ZWDGNSS

(5)

式中，ZWDGNSS為GNSS ZWD,ZWDGNSS和Π分別表示GNSS測站的天頂濕延遲與轉(zhuǎn)換系數(shù)，同式(2)一致。

由圖1(b)可知，通常情況下，MODIS像素點與GNSS測站的位置并非完全重合，因此在對MODIS PWV進行檢驗以及校正的過程中，需要確定距離GNSS站點最接近的MODIS像元，作為待檢驗的MODIS PWV。圖2為徐州地區(qū)5個測站的GNSS PWV與MODIS PWV對比的散點圖。

由圖2可以看出，GNSS PWV與MODIS PWV之間具有較好的一致性，其相關(guān)系數(shù)為0.736 4，這為利用GNSS PWV構(gòu)造MODIS PWV校正模型提供了基礎(chǔ)。根據(jù)文獻[29]所述，本文所選用的MODIS影像范圍較小，時間尺度較短，可利用GNSS PWV和MODIS PWV的實測數(shù)據(jù)，基于線性回歸分析構(gòu)造MODIS PWV的校正模型[29]，其表達式如下

PWVMODIS校正=0.406×PWVMODIS初始+31.77

(6)

式中，0.406為模型系數(shù)，31.77為常數(shù)項。表2統(tǒng)計校正前后MODIS PWV與GNSS PWV的相關(guān)系數(shù)以及精度指標(biāo)。

圖2 校正前后MODIS PWV與GNSS PWV對比Fig.2 Comparison between MODIS PWV and GNSS PWV before and after correction

由表2可知，校正后的相關(guān)系數(shù)為0.965 8。RMS從7.24 mm降低到1.35 mm，表明校正效果較好,可滿足層析試驗要求[22]。

表2 校正前后MODIS PWV與GNSS PWV的相關(guān)系數(shù)以及精度統(tǒng)計量

3 GNSS/MODIS信號緊耦合層析算法

校正后的高分辨率MODIS PWV觀測值可被添加到GNSS水汽層析模型中，但在傳統(tǒng)層析模型中通常假設(shè)一個體素塊內(nèi)水汽密度處處相等，受這一不合理假設(shè)的影響，MODIS PWV觀測值在傳統(tǒng)層析模型中只能作為觀測約束，無法充分利用MODIS觀測信號。本節(jié)基于體素節(jié)點模型對GNSS與MODIS兩種觀測信號進行參數(shù)化，提出GNSS/MODIS信號緊耦合層析算法，彌補傳統(tǒng)層析模型在融合MODIS信號方面的不足，充分發(fā)揮MODIS觀測信號的空間幾何優(yōu)勢。

3.1 融合MODIS信號的傳統(tǒng)層析模型

一般情況下，PWV觀測值可以看作沿垂直方向上的水汽密度積分值，同理，PWVMODIS可看作MODIS像素點位置處垂直方向上的水汽密度積分值，如式(7)所示

(7)

式中，ρ(s)表示垂直方向路徑上的水汽密度值，單位為g/m3，s為垂直方向路徑。在傳統(tǒng)層析模型中(圖3)，式(7)可離散化為式(8)

(8)

式中，xk表示第k層體素塊的水汽密度值，單位為g/m3；dhk表示該體素塊的高度，即層析模型第k層的厚度。

圖3 傳統(tǒng)層析模型融合MODIS信號Fig.3 Traditional tomographic model fusion with MODIS signal

在傳統(tǒng)層析模型中，受體素塊內(nèi)的水汽密度處處相等這一假設(shè)的影響，一列體素塊的水汽密度在垂直方向上的積分對應(yīng)于一個PWV觀測值。由圖1可知，MODIS PWV數(shù)據(jù)的空間分辨率為5 km×5 km，明顯高于GNSS層析模型的水平分辨率，因此會導(dǎo)致多個MODIS PWV觀測值集中分布在一列體素塊中，如圖3所示。若要將這些觀測值融合到GNSS層析模型中，需要根據(jù)一列體素塊內(nèi)的所有MODIS PWV值計算出體素塊中心位置處的PWVcenter，作為該列體素塊的PWV觀測信息?？紤]到大氣水汽分布具有一定的空間結(jié)構(gòu)性，本文采用Kriging插值來估計PWVcenter值，利用空間隨機場的結(jié)構(gòu)性信息降低估計值的不準(zhǔn)確性[30]，計算如式(9)

(9)

(10)

式中，AGNSS、AMODIS、AH以及AV分別表示GNSS信號觀測方程、MODIS信號觀測方程、水平約束方程以及垂直約束方程的系數(shù)矩陣，其中AGNSS、AMODIS表示兩種信號在體素塊層析模型中的截距信息；X為水汽密度的未知參數(shù)向量。

分析式(10)可知，在MODIS觀測條件良好情況下，MODIS信號所對應(yīng)的層析觀測方程的數(shù)量為m×n(m、n分別表示GNSS層析模型經(jīng)度方向和緯度方向上的體素塊個數(shù))。但是，由圖1可知，MODIS觀測信號極易受云朵影響，導(dǎo)致部分PWVMODIS觀測信號不可用。因此，在傳統(tǒng)體素塊層析模型中無法發(fā)揮MODIS信號空間幾何優(yōu)勢，造成有效觀測信息的浪費。

3.2 融合MODIS信號的體素節(jié)點模型

在傳統(tǒng)層析模型中，受不合理假設(shè)影響，MODIS PWV信號的利用率大幅地降低。為有效利用MODIS信號，緩解層析方程的病態(tài)性，本文對每一條MODIS信號構(gòu)造層析觀測方程，以凸顯其幾何優(yōu)勢。考慮到同一體素塊內(nèi)的多條MODIS信號到體素塊8個頂點的距離不同，采用文獻[19]提出的Trilinear parameterization方法，可以為每條MODIS PWV觀測值構(gòu)造層析觀測方程[19]。

在該方法中，將體素塊8個頂點所處位置的水汽密度設(shè)為未知參數(shù)，并用8個頂點的水汽密度的平均值表示體素塊中心的水汽密度值，以此來對傳統(tǒng)層析模型進行優(yōu)化。在本文中，考慮到大氣水汽在垂直方向分布極不均勻這一特性，垂直方向采用指數(shù)加權(quán)方法來確定中心處的水汽密度值，記為改進Trilinear parameterization方法。如圖4所示，體素塊中心O0處的水汽密度ρ(O0)可由式(11)計算出

(11)

式中，ρ(O1)和ρ(O2)分別表示上平面中心O1和下平面中心O2的水汽密度值；H為水汽標(biāo)高，取值為2 km[26]；dhk與式(8)一致。則體素中心O0的水汽密度值可用式(12)表達

ρ(P7)+ρ(P8))]

(12)

式中，ρ(Pi)表示圖4中第i個頂點處的水汽密度。

圖4 改進的trilinear parameterization方法Fig.4 The improved tilinear parameterization method

本文基于改進參數(shù)化方法分別對GNSS SWV和MODIS PWV觀測信號進行參數(shù)化。首先介紹基于體素節(jié)點模型GNSS信號的參數(shù)代過程，傳統(tǒng)GNSS層析觀測方程組的參數(shù)向量X可用新的參數(shù)向量Z代替

Xmnk,1=Mmnk,(m+1)(n+1)(k+1)·Z(m+1)(n+1)(k+1),1

(13)

式中，m、n、k分別表示傳統(tǒng)層析模型中經(jīng)度、緯度、高度3個方向上體素塊的數(shù)量；M為轉(zhuǎn)換矩陣，其詳細計算方式如式(14)所示

(14)

式中，mij為轉(zhuǎn)換矩陣M中第i行，第j列的元素，其中i=a×b×c，a、b、c分別表示傳統(tǒng)層析模型中體素塊所處的行、列和層數(shù)。基于改進Trilinear parameterization方法，GNSS觀測信號的層析方程可寫成

SWVq,1=Aq,mnk·Mmnk,(m+1)(n+1)(k+1)·

Z(m+1)(n+1)(k+1),1

(15)

式中，SWVq,1表示GNSS信號斜路徑水汽值組成的列向量。下面結(jié)合圖5介紹MODIS信號的節(jié)點參數(shù)化過程。

圖5以1條MODIS PWV觀測信號，右側(cè)對應(yīng)的是其中一個體素塊的局部放大圖。由圖1可知，該信號未穿過體素塊的中心位置，其PWV值用式(16)表達

(16)

式中，ρ(s)和s與式(7)一致，ρ(sk)表示第k層的垂直積分路徑sk上的水汽密度值。對位于第k層的分段MODIS信號而言，如圖5所示，其垂直路徑的PWV積分值可利用數(shù)值積分中的梯形公式(17)來計算[31]

(17)

式中，ρ(S0)、ρ(S1)、ρ(S2)分別表示S0、S1、S23點處的水汽密度值。因此，要構(gòu)造MODIS PWV觀測值所對應(yīng)的層析觀測方程，則需計算出S0、S1、S23點處水汽密度值與8個頂點水汽密度的數(shù)量關(guān)系。

圖5 基于體素節(jié)點模型的MODIS信號參數(shù)化過程Fig.5 The parameterization process of MODIS signal based on the voxel node tomography model

對于位于上、下層面的S1、S2來說，設(shè)其僅與同一層平面內(nèi)的4個節(jié)點的水汽密度有關(guān)，以S1為例，根據(jù)大氣水平的空間分布特性，采用高斯加權(quán)函數(shù)[3]可建立如下數(shù)量關(guān)系

ρ(S1)=φ1ρ(P1)+φ2ρ(P2)+φ3ρ(P3)+φ4ρ(P4)

(18)

對于S0而言，可以用式(11)計算其與上、下層面的S1、S2處的水汽密度的數(shù)量關(guān)系，然后利用式(18)獲得其與周圍8個頂點的空間位置關(guān)系。因此，可構(gòu)造第k層的體素塊中MODIS PWV與8個頂點的水汽密度的函數(shù)關(guān)系為

φ2ρ(P2)+φ3ρ(P3)+φ4ρ(P4))+

φ7ρ(P7)+φ8ρ(P8))]

(19)

根據(jù)式(19)可構(gòu)造一條完整的MODIS PWV所對應(yīng)的層析觀測方程，與傳統(tǒng)層析模型的融合過程相比，每一條MODIS信號可構(gòu)造一條唯一的層析觀測方程。因此，有多條MODIS信號構(gòu)成的MODIS層析觀測方程組如式(20)所示。此外，對于受云朵影響的MODIS PWV觀測值可直接舍去，并不會對MODIS層析觀測方程組造成較大的影響。

(20)

(21)

式中，BGNSS、BMODIS表示GNSS信號與MODIS信號所對應(yīng)的系數(shù)矩陣，不同于式(10)中的AGNSS和AMODIS，這兩個系數(shù)矩陣表示的是觀測信號與節(jié)點之間的空間位置關(guān)系。BH和BV分別為體素節(jié)點模型中的水平約束[3]和垂直約束[2]的系數(shù)矩陣，其構(gòu)造方法與傳統(tǒng)體素模型相似，根據(jù)節(jié)點在水平和垂直方向的位置關(guān)系分別構(gòu)造約束矩陣，Z為所有體素節(jié)點的水汽密度參數(shù)。

4 層析試驗與分析

4.1 層析方案設(shè)置

本文選取徐州地區(qū)2016年7月共15幅的MODIS影像以及相應(yīng)時段的GNSS觀測數(shù)據(jù)進行試驗。利用GAMIT/GLOBK10.6軟件對30 s采樣率的GNSS數(shù)據(jù)進行處理，ZTD和梯度項的時間分辨率為5 min。利用Saastamoinen模型[24]估計ZHD值，選用VMF1濕映射函數(shù)[25]計算SWV值。如圖1(b)所示，層析試驗區(qū)域范圍為經(jīng)度116.94°E—117.72°E，緯度33.94°N—34.64°N。在水平分辨率方面，經(jīng)度方向和緯度方向的分辨率為0.13°×0.14°；垂直分辨率方面，采用下密上疏的非等間距劃分方法[13]，共劃分為15層。并且根據(jù)徐州地區(qū)10年(2009—2018年)每年7月探空數(shù)據(jù)提供的水汽廓線，確定最優(yōu)對流層頂高度為11 km，考慮到不同天氣條件的影響，圖6展示了2009—2018年徐州地區(qū)每年7月的水汽垂直分布特征。此外，本文采用代數(shù)重構(gòu)算法解算層析觀測方程組，利用層析時段前3天的探空垂直廓線信息作為迭代初值進行層析試驗[4]。

圖6 徐州地區(qū)2009—2018年每年7月水汽垂直分布圖Fig.6 Vertical distribution of water vapor in Xuzhou during each July from 2009 to 2018

為了驗證本文提出的緊耦合算法，共設(shè)計3種試驗方案來進行對比分析。

方案1：采用傳統(tǒng)層析模型，僅利用GNSS數(shù)據(jù)作為觀測信息。

方案2：采用傳統(tǒng)層析模型，并利用GNSS和MODIS數(shù)據(jù)構(gòu)造層析觀測方程進行解算。

方案3：采用體素節(jié)點模型，利用改進Trilinear parameterization方法分別對GNSS信號與MODIS信號進行參數(shù)化，并利用兩者構(gòu)成的層析觀測方程組進行層析解算。

4.2 MODIS PWV信號貢獻性分析

在體素節(jié)點模型中加入MODIS信號后，部分未被GNSS信號線穿刺的體素塊將會被MODIS信號線穿過，可以有效彌補傳統(tǒng)層析模型幾何結(jié)構(gòu)的缺陷。圖7展示了加入MODIS信號前后，層析模型中第1層、第5層、第10層、第15層的體素塊所穿過信號的平均數(shù)量對比。

圖7 加入MODIS信號前后穿過體素塊的信號數(shù)量對比Fig.7 The comparison of the number of observation signals passing through the voxels before and after adding MODIS signals

由圖7可知，一方面，在0～2 km的近地層，傳統(tǒng)層析模型中有大量的體素塊未被GNSS信號穿過，尤其第1層中只有GNSS測站所在的體素塊被觀測信號穿刺，加入MODIS信號后，這些空白體素塊均有信號線通過，較好地改善了層析模型的觀測幾何結(jié)構(gòu)。另一方面，在4 km以上的高層中，加入MODIS信號之后同樣使得每個體素塊所穿過的信號數(shù)量增加4～6條，為層析模型引入更多的有用觀測信息。圖8統(tǒng)計了每個層析時段內(nèi)有效的GNSS信號和MODIS信號的數(shù)量，將MODIS信號添加到層析系統(tǒng)后，平均有效觀測信號數(shù)量從366條增加到491條，提高了34.15%，這表明引入MODIS數(shù)據(jù)可以提高有效觀測值數(shù)量。

圖8 層析時段內(nèi)不同類型觀測信號的數(shù)量對比Fig.8 Comparison of the number of different observation signals during the tomography period

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

在層析結(jié)果精度評定中，通常以探空站在UTC 00:00和UTC 12:00兩個時刻提供的高精度垂直水汽廓線信息為參考值，但本文中的層析時段無法與探空觀測時間完全一致，因此選用距層析時段最近的探空數(shù)據(jù)對水汽反演結(jié)果進行檢驗[21-22]。圖9(a)和圖9(b)分別給出了層析時段內(nèi)3種方案的層析結(jié)果的均方根誤差(root mean squared error，RMSE)和平均絕對誤差(mean square error,MAE)的對比情況。由圖9可知，在大多數(shù)層析歷元內(nèi)，方案3的RMSE和MAE值都要優(yōu)于方案1和方案2，方案2的統(tǒng)計量略優(yōu)于方案1。此外，表3列出了15個層析時段內(nèi)不同方案的層析結(jié)果與探空數(shù)據(jù)對比的統(tǒng)計結(jié)果，可知3種方案的層析結(jié)果的平均RMSE分別為2.43、2.24、1.82 g/m3；平均MAE分別為1.96、1.72、1.43 g/m3。相對于傳統(tǒng)層析模型(方案1)，本文提出的算法的反演精度(RMS和MAE)平均提高了26.07%，表明融合MODIS信號的體素節(jié)點模型具有更好的反演能力。

圖9 層析時段內(nèi)3種方案的層析結(jié)果的RMSE和MAE對比Fig.9 Comparison of the RMSE and MAE of the tomographic results derived from three schemes during the tomography periods

表3 15個層析時段內(nèi)不同方案的層析結(jié)果與探空數(shù)據(jù)對比的統(tǒng)計結(jié)果

為了進一步比較不同算法反演水汽的時空分布情況，對兩種天氣條件下(晴天和雨天)反演的水汽廓線進行對比分析。其中，7-01、7-11、7-21、7-25為晴天；7-17、7-26、7-27、7-28對應(yīng)的是雨天。圖11給出了兩種天氣條件下3種方案反演得到的水汽廓線與探空站提供的水汽廓線的對比圖。對比來看，方案1和方案2得到的水汽廓線較為接近，說明將MODIS信號融合到體素塊層析模型中并不能起到有效作用，但方案3的水汽廓線明顯優(yōu)于兩者，并與探空站反演的水汽廓線相一致。

此外，根據(jù)水汽廓線比對情況，本節(jié)對層析時段內(nèi)不同高度層的水汽密度進行精度評定，圖11給出了3種方案的層析結(jié)果與探空站所在列體素塊的水汽密度的誤差對比圖。可以看出，在4 km以上的高度層中，優(yōu)化算法的層析結(jié)果質(zhì)量優(yōu)于傳統(tǒng)算法；在2～4 km的高度層中，方案3的反演精度略優(yōu)于前兩種方案；0～2 km的近地層中，富含60%的大氣水汽含量，是造成信號延遲的主要區(qū)域。通過定量計算可得3種算法在0～2 km的層析結(jié)果平均RMSE分別為3.87 g/m3、3.75 g/m3、2.55 g/m3，可以看出優(yōu)化算法的反演精度明顯高于傳統(tǒng)算法，進一步表明融合MODIS觀測信號的可改善近地層三維水汽場的重構(gòu)質(zhì)量。

圖10 3種方案在不同天氣條件下的層析水汽廓線對比Fig.10 Comparison of tomographic water vapor profiles derived from three schemes under different weather conditions

圖11 層析時段內(nèi)3種方案層析結(jié)果與探空參考值的誤差對比Fig.11 The difference between tomographic results and radiosonde along the radiosonde vertical columns

5 結(jié) 論

本文基于高分辨率MODIS PWV數(shù)據(jù)，提出GNSS/MODIS信號緊耦合水汽層析算法。憑借體素節(jié)點模型的建模優(yōu)勢，分別建立基于GNSS信號和MODIS信號的層析觀測方程組。利用2016年7月徐州地區(qū)15幅MODIS影像和同一時刻下的GNSS實測數(shù)據(jù)，系統(tǒng)地比較分析了緊耦合算法與傳統(tǒng)算法的層析結(jié)果的質(zhì)量及精度。主要結(jié)論如下：

(1) 引入MODIS觀測信號后，平均有效觀測信號數(shù)量提高了34.15%，且所有體素塊中均有觀測信號穿過，可有效改善傳統(tǒng)層析模型的GNSS信號幾何結(jié)構(gòu)缺陷。

(2) 與傳統(tǒng)算法的層析結(jié)果相比，融合MODIS信號的傳統(tǒng)算法和所提算法層析結(jié)果的平均RMSE值由2.43 g/m3分別降低為2.24 g/m3和1.82 g/m3，平均反演精度分別提高了7.81%和25.10%；與融合MODIS信號的傳統(tǒng)算法相比，緊耦合算法的反演精度提高18.75%。

(3) 在0～2 km的近地層，緊耦合算法反演的三維水汽場質(zhì)量明顯優(yōu)于其他地方兩種層析結(jié)果，該高度層的平均RMSE值由3.87 g/m3(方案1)和3.75 g/m3(方案2)降低為2.55 g/m3(方案3)，反演精度分別提高了34.11%和32%，這表明緊耦合算法可最大化利用MODIS觀測信號，有效改善層析結(jié)果質(zhì)量。

致謝：特別感謝徐州市自然資源和規(guī)劃局提供的XZCORS網(wǎng)監(jiān)測數(shù)據(jù)，美國國家航空航天局提供的MODIS影像，美國懷俄明大學(xué)提供的長期探空數(shù)據(jù)，麻省理工學(xué)院提供的GAMIT/GLOBK軟件。