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

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

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

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

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

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

1 GNSS水汽層析原理

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

(1)

SWV=Π·SWD

(2)

(3)

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

(4)

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

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

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

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

MODIS數(shù)據(jù)包括Terra衛(wèi)星平臺(tái)提供的MOD05_L2以及Aqua衛(wèi)星提供的MYD05_L2兩種水汽產(chǎn)品，包含每個(gè)像元的經(jīng)緯度、掃描時(shí)間、太陽天頂角和方位角、傳感器天頂角和方位角、近紅外水汽值、紅外水汽值、云掩膜產(chǎn)品等13個(gè)參數(shù)[24]。本文采用來自Terra平臺(tái)的MOD05_L2水汽產(chǎn)品進(jìn)行層析試驗(yà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ù)會(huì)成倍地增加GNSS層析模型的計(jì)算量，因此選用5 km×5 km分辨率的水汽產(chǎn)品用于層析試驗(yàn)研究。圖1(a)和(b)分別顯示了2016年7月1日14:00時(shí)我國(guó)東部地區(qū)的MODIS水汽數(shù)據(jù)以及徐州地區(qū)的局部放大圖，圖1(a)中的灰色區(qū)域表示被云朵污染的部分。

2.2 MODIS PWV檢驗(yàn)與校正

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

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

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

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

PWV=Π·ZWDGNSS

(5)

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

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

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

PWVMODIS校正=0.406×PWVMODIS初始+31.77

(6)

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

圖2 校正前后MODIS PWV與GNSS PWV對(duì)比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，表明校正效果較好,可滿足層析試驗(yàn)要求[22]。

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

3 GNSS/MODIS信號(hào)緊耦合層析算法

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

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

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

(7)

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

(8)

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

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

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

(9)

(10)

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

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

3.2 融合MODIS信號(hào)的體素節(jié)點(diǎn)模型

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

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

(11)

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

ρ(P7)+ρ(P8))]

(12)

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

圖4 改進(jìn)的trilinear parameterization方法Fig.4 The improved tilinear parameterization method

本文基于改進(jìn)參數(shù)化方法分別對(duì)GNSS SWV和MODIS PWV觀測(cè)信號(hào)進(jìn)行參數(shù)化。首先介紹基于體素節(jié)點(diǎn)模型GNSS信號(hào)的參數(shù)代過程，傳統(tǒng)GNSS層析觀測(cè)方程組的參數(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個(gè)方向上體素塊的數(shù)量；M為轉(zhuǎn)換矩陣，其詳細(xì)計(jì)算方式如式(14)所示

(14)

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

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

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

(15)

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

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

(16)

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

(17)

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

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

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

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

(18)

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

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

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

(19)

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

(20)

(21)

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

4 層析試驗(yàn)與分析

4.1 層析方案設(shè)置

本文選取徐州地區(qū)2016年7月共15幅的MODIS影像以及相應(yīng)時(shí)段的GNSS觀測(cè)數(shù)據(jù)進(jìn)行試驗(yàn)。利用GAMIT/GLOBK10.6軟件對(duì)30 s采樣率的GNSS數(shù)據(jù)進(jìn)行處理，ZTD和梯度項(xiàng)的時(shí)間分辨率為5 min。利用Saastamoinen模型[24]估計(jì)ZHD值，選用VMF1濕映射函數(shù)[25]計(jì)算SWV值。如圖1(b)所示，層析試驗(yàn)區(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)對(duì)流層頂高度為11 km，考慮到不同天氣條件的影響，圖6展示了2009—2018年徐州地區(qū)每年7月的水汽垂直分布特征。此外，本文采用代數(shù)重構(gòu)算法解算層析觀測(cè)方程組，利用層析時(shí)段前3天的探空垂直廓線信息作為迭代初值進(jìn)行層析試驗(yàn)[4]。

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

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

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

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

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

4.2 MODIS PWV信號(hào)貢獻(xiàn)性分析

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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