王 昊，丁 楠，張文淵，馮遵德，趙長(zhǎng)勝，閆香蓉

1.江蘇師范大學(xué)地理測(cè)繪與城鄉(xiāng)規(guī)劃學(xué)院，江蘇 徐州 221116；2.中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院，江蘇 徐州 221116

大氣水汽是影響氣候變化的關(guān)鍵氣象要素，對(duì)大氣能量輸送和對(duì)流天氣的演變具有重要影響作用[1-2]。傳統(tǒng)水汽探測(cè)手段，如無線電探空技術(shù)、水汽微波輻射計(jì)等，由于具有諸多局限性，已不能滿足目前實(shí)時(shí)、快捷、精確地獲取水汽要素的要求。隨著GNSS技術(shù)在氣象學(xué)中的應(yīng)用發(fā)展，利用GNSS技術(shù)進(jìn)行對(duì)流層三維水汽場(chǎng)的反演已快速成為一種新興水汽探測(cè)手段。GNSS水汽層析技術(shù)具有精確、快速、全天候和低成本等諸多優(yōu)勢(shì)[3]，已在大氣水汽監(jiān)測(cè)中得到越來越多的應(yīng)用。

文獻(xiàn)[4]率先發(fā)展了GNSS水汽層析技術(shù)，研究了垂直分層對(duì)水汽層析結(jié)果精度的影響規(guī)律，基于夏威夷Kilauea GPS網(wǎng)進(jìn)行了水汽層析試驗(yàn)，通過對(duì)比歐洲中期天氣預(yù)報(bào)中心(European Centre for Medium-Range Weather Forecasts，ECMWF)的數(shù)據(jù)，得出了垂直分層的最小分辨率需大于300 m以削弱噪聲影響并保證層析結(jié)果精度的結(jié)論。垂直分層是GNSS水汽層析的關(guān)鍵技術(shù)之一，基于文獻(xiàn)[4]的水汽層析思路，諸多學(xué)者對(duì)這一關(guān)鍵技術(shù)進(jìn)行了研究，目前主要分為兩類：一類為采用恒定高度間隔的均勻分層；另一類為采用不同高度間隔的非均勻分層。在均勻分層的研究中，文獻(xiàn)[5]采用3種不同高度間隔(600、800、1000 m)的等距均勻垂直分層進(jìn)行層析試驗(yàn)，發(fā)現(xiàn)采用均勻分層時(shí)以800 m為間隔最優(yōu)。在非均勻分層的研究中，文獻(xiàn)[6]以四川省為例對(duì)非均勻分層與均勻分層進(jìn)行對(duì)比試驗(yàn)，研究并分析了采用非均勻分層對(duì)層析結(jié)果的影響。文獻(xiàn)[7]基于探空數(shù)據(jù)研究了層析區(qū)域的頂層邊界高度，提出了兩種不同非均勻分層模型，并對(duì)比了兩種模型的衛(wèi)星信號(hào)利用率和體素塊利用率。文獻(xiàn)[8]提出了一種體素塊優(yōu)化的層析區(qū)域建模思路，在垂直方向上以1973—2012年探空水汽廓線為依據(jù)對(duì)香港地區(qū)進(jìn)行4種不同的垂直分層建模，對(duì)比分析了不同高度間隔的非均勻分層模型對(duì)層析結(jié)果質(zhì)量的影響。

傳統(tǒng)的垂直分層建模方法將研究區(qū)域按等高度間距于垂直方向進(jìn)行劃分[9]，該方法屬于均勻分層，優(yōu)勢(shì)在于操作簡(jiǎn)便且分層后獲得的離散區(qū)域形狀規(guī)則[10]。由于大氣水汽在垂直方向上呈指數(shù)遞減趨勢(shì)，因此建模時(shí)采用均勻分層的不足之處在于均勻的垂直分辨率不符合大氣水汽的實(shí)際分布，會(huì)引入一定的誤差。而非均勻分層模型的垂直分辨率更加符合大氣水汽的實(shí)際分布，采用非均勻分層模型進(jìn)行離散化能有效減小上述誤差。為削弱垂直分層建模過程中誤差的影響，非均勻分層模型已被廣泛運(yùn)用于GNSS水汽層析技術(shù)之中[11-12]，但目前非均勻分層的建模方法仍然存在不足，非均勻分層模型的垂直分辨率設(shè)置尚未形成科學(xué)可靠且適用性強(qiáng)的參考標(biāo)準(zhǔn)，部分層析試驗(yàn)在采用非均勻分層模型時(shí)僅根據(jù)固定的經(jīng)驗(yàn)設(shè)置上疏下密的分層高度區(qū)間，并無科學(xué)可靠的分層依據(jù)[13]；另一部分試驗(yàn)雖基于先驗(yàn)水汽密度進(jìn)行建模，但仍存在分層模型適用性不高的問題，在層析區(qū)域變更或大氣水汽垂直分布特征發(fā)生變化的情況下無法實(shí)現(xiàn)方便快捷的自適應(yīng)分層[7-8,14]。

針對(duì)上述問題，本文提出一種自適應(yīng)非均勻指數(shù)分層(adaptive non-uniform exponential vertical stratification，ANEVS)的GNSS水汽層析方法，詳細(xì)闡述了該垂直分層方法的分層步驟和具體公式。該方法基于層析區(qū)域內(nèi)大氣水汽于垂直方向上呈指數(shù)遞減的分布特征，能夠自適應(yīng)地完成任意給定層析區(qū)域的非均勻垂直分層建模，優(yōu)化各層水汽密度參數(shù)的分配，縮小底層各層水汽密度參數(shù)之間的差異，提高了垂直分層的合理性和真實(shí)性，進(jìn)而改善層析結(jié)果質(zhì)量，提高層析解算精度，同時(shí)彌補(bǔ)了垂直分層模型適用性低的問題。

1 GNSS水汽層析基本原理

GNSS水汽層析技術(shù)將連續(xù)的研究區(qū)域離散化為數(shù)量有限的體素塊，以衛(wèi)星信號(hào)的斜路徑濕延遲(slant wet delay，SWD)為觀測(cè)值構(gòu)建大氣水汽參數(shù)、信號(hào)傳播距離與信號(hào)傳播路徑上的總水汽含量[15]之間的函數(shù)關(guān)系，反演出水汽的三維空間分布信息[16]。

GNSS水汽層析以每條信號(hào)的SWD為觀測(cè)值。忽略信號(hào)傳播過程中的彎曲量，假設(shè)傳播路徑為三維空間直線，同時(shí)假設(shè)各體素塊內(nèi)水汽密度在一個(gè)層析歷元內(nèi)為常數(shù)。在此基礎(chǔ)上建立模型為

(1)

式中，SWDi表示第i條衛(wèi)星信號(hào)的斜路徑濕延遲；lij表示第i條信號(hào)穿過第j個(gè)體素塊的直線長(zhǎng)度；Dj表示第j個(gè)體素塊內(nèi)的水汽參數(shù)，這里為該體素塊的待求水汽密度；n表示總體素塊個(gè)數(shù)。

基于上述模型構(gòu)建的觀測(cè)方程組系數(shù)矩陣嚴(yán)重秩虧，因此在構(gòu)建層析方程組時(shí)，一般通過加入約束條件以解決層析方程組的不適定性問題[17-18]。約束條件主要包括垂直約束、水平約束和頂層約束。垂直約束采用指數(shù)函數(shù)進(jìn)行約束，其模型為[19]

(2)

式中，hi和hj分別表示大氣高度；ρi和ρj分別表示高度hi和hj處的水汽密度；H表示水汽標(biāo)高[20]。

水平約束采用高斯加權(quán)函數(shù)進(jìn)行約束，其模型為[21]

w1D1+w2D2+…+wi-1Di-1+wi+1Di+1+

+wjDj=Di

(3)

(4)

式中，Di表示同層第i個(gè)體素塊內(nèi)的水汽密度；wi表示同層第i個(gè)水汽密度的權(quán)系數(shù)；dij表示第i個(gè)體素塊與同層第j個(gè)體素塊中點(diǎn)之間的距離；σ表示平滑因子[22]，其取值與體素塊水平方向的長(zhǎng)度有關(guān)[21]，本文的平滑因子取值為15。

頂層約束將最頂層體素塊內(nèi)的水汽密度按照經(jīng)驗(yàn)強(qiáng)制約束為某一定值，部分層析試驗(yàn)將頂層水汽密度值約束為0.01 g/m3[23]，但頂層體素塊實(shí)際具有一定高度，僅頂層體素塊頂端水汽密度值接近0.01 g/m3，中間及底端水汽密度大于0.01 g/m3，因此本文以先驗(yàn)探空數(shù)據(jù)為依據(jù)，采用頂層體素塊高度區(qū)間內(nèi)的探空水汽密度均值為約束值。

結(jié)合觀測(cè)方程和約束條件的層析方程組可以表示為

(5)

式中，A0、AV、AH、AT分別表示觀測(cè)值系數(shù)矩陣、垂直約束系數(shù)矩陣、水平約束系數(shù)矩陣、頂層約束系數(shù)矩陣；DW表示待求水汽密度向量；b0、bV、bH、bT分別表示對(duì)應(yīng)的觀測(cè)值向量，其中，bV、bH一般情況下為0向量；方程總個(gè)數(shù)m等于各部分方程個(gè)數(shù)o、t、v、h之和。

式(5)的解算方法有奇異值分解法、最小二乘法、代數(shù)重構(gòu)法(algebraic reconstruction technique，ART)等[24-25]，本文采用ART對(duì)層析方程組進(jìn)行解算。ART避免了系數(shù)矩陣的求逆運(yùn)算，具有抗噪性強(qiáng)的優(yōu)點(diǎn)，即使在觀測(cè)方程較少、信號(hào)分布不均勻和信號(hào)噪聲較大的不利條件下，也能解算出精度較高的結(jié)果[26-27]。

2 基于水汽垂直分布特征的自適應(yīng)非均勻指數(shù)分層方法

自適應(yīng)非均勻指數(shù)分層方法總體分為兩大步驟：首先，對(duì)層析區(qū)域進(jìn)行大氣水汽的垂直分布建模，以該區(qū)域的先驗(yàn)水汽廓線信息為依據(jù)獲得大氣水汽的垂直分布特征函數(shù)；然后，基于該水汽垂直分布特征函數(shù)，通過迭代的方式進(jìn)行各層高度區(qū)間的計(jì)算，實(shí)現(xiàn)對(duì)任意層析區(qū)域進(jìn)行任意層數(shù)劃分的自適應(yīng)分層。

2.1 水汽垂直分布建模

非均勻指數(shù)分層基于先驗(yàn)水汽垂直分布特征函數(shù)，以水汽密度區(qū)間計(jì)算分層高度區(qū)間，在進(jìn)行自適應(yīng)分層前需先對(duì)層析區(qū)域大氣水汽的垂直分布進(jìn)行建模，以獲取高度與水汽密度之間的函數(shù)關(guān)系。

基于水汽密度隨高度增加呈指數(shù)遞減規(guī)律對(duì)大氣水汽的垂直分布特征進(jìn)行建模[28]，高度h處的水汽密度表達(dá)為

(6)

式中，ρh表示高度h處的水汽密度，單位為g/m3；ρ0表示地表處水汽密度參數(shù)；Dv表示水汽垂直分布特征參數(shù)；H表示水汽標(biāo)高，單位為m[20]。

由式(6)可知，某一高度h處的水汽密度值ρh(g/m3)可表示為高度h的函數(shù)，其中包含水汽密度參數(shù)ρ0和水汽垂直分布特征參數(shù)Dv兩個(gè)未知參數(shù)，通過對(duì)研究區(qū)域的先驗(yàn)水汽密度數(shù)據(jù)進(jìn)行擬合可獲得式(6)中的兩個(gè)未知參數(shù)。獲得先驗(yàn)水汽密度數(shù)據(jù)的方法有多種，以探空氣象站的探空數(shù)據(jù)為例，探空數(shù)據(jù)具有高垂直分辨率的優(yōu)點(diǎn)，高精度的先驗(yàn)探空水汽密度可通過原始探空資料中包含的溫度、露點(diǎn)溫度等觀測(cè)數(shù)據(jù)進(jìn)行計(jì)算，水汽密度Dw計(jì)算公式為

(7)

式中，e表示水汽分壓，單位為hPa；Rv=461.5 J·kg-1·K-1；T表示溫度，單位為K。其中，e可由世界氣象組織推薦的Goff-Gratch公式計(jì)算

(8)

式中，T1=273.16 K；Td表示露點(diǎn)溫度，單位為K；e表示水汽分壓，單位為hPa。

上述部分詳細(xì)闡述了利用探空數(shù)據(jù)計(jì)算先驗(yàn)水汽密度的公式和步驟，但由于探空數(shù)據(jù)的空間分辨率和時(shí)間分辨率較低，可能出現(xiàn)層析區(qū)域或?qū)游鰰r(shí)間段內(nèi)探空數(shù)據(jù)無法獲取的情況。針對(duì)該情況，利用ECMWF或美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心(European Centre for Medium-Range Weather Forecasts，NCEP)提供的再分析資料中溫度和相對(duì)濕度等氣象數(shù)據(jù)同樣可以獲得可靠的先驗(yàn)水汽密度數(shù)據(jù)[18,29]。

圖1為水汽垂直分布特征擬合圖，其中散點(diǎn)代表每一高度位置處所對(duì)應(yīng)的實(shí)測(cè)探空水汽密度，曲線代表最小二乘擬合后的函數(shù)曲線，設(shè)擬合后曲線的函數(shù)表達(dá)式為

圖1 水汽垂直分布特征擬合Fig.1 Vertical distribution characteristics of water vapor

(9)

式中，x表示高度，單位為m；y表示水汽密度，單位為g/m3；a和b分別表示地表處水汽密度參數(shù)ρ0和水汽垂直分布特征參數(shù)Dv擬合后的具體數(shù)值；H表示水汽標(biāo)高，單位為m。式(9)反映了層析區(qū)域內(nèi)大氣水汽于垂直方向上的總體分布特征，為下一步自適應(yīng)分層提供了高度值與水汽密度值之間的具體函數(shù)關(guān)系。

2.2 自適應(yīng)非均勻指數(shù)分層

在獲得如式(9)所示的層析區(qū)域水汽垂直分布特征函數(shù)后，可通過如下方法進(jìn)行自適應(yīng)垂直分層：將層析區(qū)域劃分為L(zhǎng)層，基于文獻(xiàn)[4]的研究，垂直分層間隔需大于300 m[4]，首先將第1層高度區(qū)間設(shè)置為[0,300]，根據(jù)式(9)計(jì)算300 m處水汽密度為a·eb·(-300/H)(g/m3)，第1層水汽密度區(qū)間為[a,a·eb·(-300/H)]，現(xiàn)還需將剩余高度劃分為L(zhǎng)-1個(gè)高度區(qū)間，計(jì)算剩余L-1層每層水汽密度間隔為(a·eb·(-300/H)-0.1)/(L-1)，單位為m，以該值計(jì)算各層預(yù)設(shè)水汽密度區(qū)間，得剩余L-1層的預(yù)設(shè)水汽密度區(qū)間表達(dá)式為

(10)

式中，a、b、H表示含義與式(9)相同；n表示層數(shù)(2≤n≤L)；L表示預(yù)計(jì)劃分總層次數(shù)。

將剩余L-1層對(duì)應(yīng)式(10)的預(yù)設(shè)水汽密度區(qū)間代入式(9)，反推計(jì)算出剩余L-1層的預(yù)設(shè)高度區(qū)間，其中第2層預(yù)設(shè)高度區(qū)間為

(11)

式中，a、b、H、L表示含義與式(10)相同。

若該層高度間隔大于300 m，則以此時(shí)的預(yù)設(shè)高度區(qū)間為結(jié)果，完成非均勻指數(shù)分層；若該層高度間隔小于300 m，則將該層高度區(qū)間繼續(xù)以最小高度間隔300 m設(shè)置為1層，例如第2層高度區(qū)間為[300,600]，密度區(qū)間為[a·eb·(-300/H),a·eb·(-600/H)]，剩余層數(shù)為L(zhǎng)-2，每層密度間隔為(a·eb·(-600/H)-0.1)/(L-2)，以新的密度間隔計(jì)算新的剩余層預(yù)設(shè)密度區(qū)間，在進(jìn)行第m次以300 m為高度間隔劃分最底層后，剩余層的預(yù)設(shè)密度區(qū)間可表示為

(12)

式中，a、b、H表示含義與式(11)相同；m表示以300 m為高度間隔進(jìn)行劃分的層次數(shù)量；n表示剩余第m+1層至L層的層號(hào)(m+1≤n≤L)；L表示預(yù)計(jì)劃分總層次數(shù)。

每次迭代，將式(12)的剩余層預(yù)設(shè)密度區(qū)間代入式(9)反推新的預(yù)設(shè)高度區(qū)間，判斷其中最低層高度間隔是否大于300 m，最低層高度間隔表示為

(13)

式中，a、b、H、m、L含義與式(12)相同；Hi表示剩余L-m層中最底層高度間隔，單位為m，若Hi大于300 m則完成非均勻指數(shù)分層，反之則繼續(xù)以300 m為高度間隔劃分剩余L-m層中的最底層并以此類推。

上述的自適應(yīng)非均勻指數(shù)分層如圖2所示。ANEVS進(jìn)行垂直分層時(shí)，先通過所需分層層數(shù)計(jì)算水汽密度區(qū)間，再通過水汽密度區(qū)間和先驗(yàn)水汽垂直分布特征函數(shù)反推高度區(qū)間，通過迭代判斷并調(diào)整最底層高度間隔。該方法能根據(jù)層析區(qū)域的大氣水汽垂直分布特征和所需分層層數(shù)自動(dòng)調(diào)整分層高度區(qū)間，同時(shí)自適應(yīng)地控制各層水汽密度差異在一個(gè)很小的范圍內(nèi)，做到細(xì)分底層區(qū)域密度不同的水汽、合并高層區(qū)域密度相近的水汽，保證了垂直分層的精度和真實(shí)性，且該方法適用于任意層析區(qū)域，即使在水汽垂直分布特征發(fā)生改變的情況下，也能迅速根據(jù)改變后的水汽垂直分布特征自動(dòng)進(jìn)行所需分層數(shù)的自適應(yīng)非均勻分層，具有極強(qiáng)的適用性。

圖2 自適應(yīng)非均勻指數(shù)分層 Fig.2 Adaptive non-uniform exponential vertical stratification

3 GNSS水汽層析試驗(yàn)驗(yàn)證

3.1 試驗(yàn)區(qū)域及時(shí)段

本文以香港為試驗(yàn)區(qū)域，圖3為試驗(yàn)區(qū)域內(nèi)參考站及探空氣象站分布圖，其中實(shí)心圓圖形表示19個(gè)參考站所在位置，實(shí)心三角形圖形表示探空氣象站所在位置。本文利用19個(gè)連續(xù)運(yùn)行參考站(continuously operating reference system，CORS)的GPS觀測(cè)數(shù)據(jù)進(jìn)行水汽層析試驗(yàn)，以京士柏探空氣象站(HKKP)的探空(radiosonde，RS)數(shù)據(jù)為參考，將HKKP位置天頂方向的層析水汽密度與探空水汽密度進(jìn)行對(duì)比分析。

圖3 香港地區(qū)參考站及探空站地理位置Fig.3 Geographic distribution of the Hong Kong SatRef with GNSS reference stations and radiosonde

試驗(yàn)時(shí)段選取2019年8月(31d)，試驗(yàn)時(shí)段內(nèi)香港地區(qū)正值夏季，氣候炎熱，且由于受熱帶氣旋“韋帕”及“白鹿”的影響，共發(fā)生3次大規(guī)模降雨，因此香港2019年8月同時(shí)具備了暴雨和炎熱無雨的不同水汽情況，十分適宜進(jìn)行水汽層析試驗(yàn)。

3.2 試驗(yàn)方案設(shè)置

觀測(cè)值SWD的獲取使用GAMIT10.7對(duì)GPS觀測(cè)數(shù)據(jù)進(jìn)行處理，為消除多路徑效應(yīng)的影響，將衛(wèi)星截止高度角設(shè)置為15°；引入固體潮模型和海潮模型進(jìn)行改正；映射函數(shù)模型設(shè)置為VMF1；引入BJFS、CHAN、USUD 3個(gè)外網(wǎng)輔助IGS站進(jìn)行聯(lián)合解算以削弱天頂對(duì)流層延遲(zenith total delay，ZTD)的相關(guān)性。

層析方程組的解算采用ART技術(shù)。該技術(shù)中迭代初值的精度對(duì)層析結(jié)果的準(zhǔn)確性具有重要影響，高精度的迭代初值可同時(shí)提高迭代速度和層析結(jié)果質(zhì)量。本文根據(jù)香港地區(qū)2010—2018年探空數(shù)據(jù)計(jì)算獲得的不同年份8月全部垂直水汽廓線均值作為迭代初值進(jìn)行層析試驗(yàn)[27]。

由于在試驗(yàn)時(shí)段內(nèi)只能獲取每日UTC 00∶00時(shí)和UTC 12∶00時(shí)的探空數(shù)據(jù)作為參考，故每日只選取UTC 00∶00時(shí)和UTC 12∶00時(shí)兩個(gè)時(shí)段進(jìn)行水汽層析試驗(yàn)。本文以水汽密度等于0.1 g/m3的高度作為對(duì)流層層頂，通過試驗(yàn)時(shí)段內(nèi)62次探空頂層高度的平均值確定試驗(yàn)區(qū)域頂層高度，經(jīng)計(jì)算頂層高度平均值為10 770 m，標(biāo)準(zhǔn)差(standard deviation，STD)為935 m。水汽層析試驗(yàn)區(qū)域水平方向范圍為22.204°N—22.564°N與113.844°E—114.384°E之間，水平分辨率為0.090°×0.090°[23]，垂直方向范圍為地面至10 770 m之間，垂直分辨率由分層方案決定。

方案A：該方案基于傳統(tǒng)均勻分層，均勻分層高度間隔在800 m時(shí)，可以獲得較好的層析結(jié)果[5]，因此該方案將10 770 m高的層析區(qū)域分為13層，使每層間隔為828.5 m，最接近800 m。頂層約束的先驗(yàn)水汽密度以先驗(yàn)探空數(shù)據(jù)為依據(jù)，由于頂層約束值作為層析方程組的觀測(cè)值直接參與三維水汽密度場(chǎng)的反演，不能采用試驗(yàn)時(shí)段內(nèi)的探空數(shù)進(jìn)行先驗(yàn)水汽密度的計(jì)算，且大氣水汽存在周期性變化規(guī)律[30-33]，故采用試驗(yàn)?zāi)攴萸耙荒?2018年)同時(shí)段(夏季8月)的探空數(shù)據(jù)計(jì)算先驗(yàn)水汽密度，經(jīng)計(jì)算，水汽密度均值為0.23 g/m3，故將頂層約束設(shè)置為0.23 g/m3。

方案B：該方案基于ANEVS，分層數(shù)設(shè)置為與方案A相同的13層。大氣水汽的垂直分布特征函數(shù)采用先驗(yàn)探空水汽密度數(shù)據(jù)進(jìn)行擬合。由于探空氣象數(shù)據(jù)具有高精度、高垂直分辨率的優(yōu)點(diǎn)，同時(shí)，以探空站點(diǎn)處的水汽密度信息代表探空站周圍一定區(qū)域內(nèi)大氣水汽密度真值的可靠性已經(jīng)被證實(shí)[34-36]，故探空站的垂直水汽密度廓線常被用作層析區(qū)域內(nèi)水汽密度的參考真值，以輔助分析水汽層析的反演結(jié)果[7-8,37]。本文的層析試驗(yàn)區(qū)域?qū)儆谛》秶鷮游鰠^(qū)域，基于高精度的先驗(yàn)探空水汽密度建立的水汽垂直分布模型在該層析區(qū)域內(nèi)具備有效性。由于水汽垂直分布特征函數(shù)直接影響了觀測(cè)方程系數(shù)矩陣的數(shù)值，因此作為先驗(yàn)信息的水汽垂直分布特征函數(shù)不能采用試驗(yàn)時(shí)段內(nèi)的探空水汽密度數(shù)據(jù)進(jìn)行擬合，基于大氣水汽的周期性變化特征[30-33]，以試驗(yàn)?zāi)攴萸耙荒?2018年)同時(shí)段(夏季8月)的探空水汽密度作為先驗(yàn)信息在時(shí)間域上具備代表性，擬合后的指數(shù)函數(shù)表達(dá)式為

(14)

式中，Dw表示水汽密度值，單位為g/m3；h表示高度，單位為m。探空數(shù)據(jù)擬合后的R2為0.977 4；均方根誤差(root mean squared error，RMSE)為0.919 g/m3?；谑?14)的水汽垂直分布特征函數(shù)對(duì)試驗(yàn)區(qū)域運(yùn)用ANEVS進(jìn)行垂直分層后，將0～1500 m以300 m為高度間隔劃分5層，1500～10 770 m劃分為剩余8層，每層水汽密度間隔約為1.70 g/m3，第6層為1500～1838 m，高度間隔338 m，大于300 m，經(jīng)前一年同時(shí)段探空數(shù)據(jù)計(jì)算頂層水汽密度均值為0.87 g/m3，故將頂層約束設(shè)置為0.87 g/m3。

分別以方案A和方案B對(duì)試驗(yàn)區(qū)域進(jìn)行垂直分層建模，圖4為兩種分層方案的分層建模示意圖。

圖4 方案A及方案B垂直分層模型的垂直分辨率設(shè)置Fig.4 The vertical resolution of the vertical stratification model of project A and project B

4 GNSS水汽層析試驗(yàn)結(jié)果分析

4.1 不同垂直分層方案的解算精度分析

圖5為試驗(yàn)時(shí)段內(nèi)UTC 00∶00時(shí)與UTC 12∶00時(shí)層析歷元下兩種方案的單天RMSE值對(duì)比圖。其中方案A的RMSE數(shù)值最大出現(xiàn)在DOY 213 UTC 00∶00時(shí)(2.15 g/m3)，數(shù)值最小出現(xiàn)在DOY 214 UTC 12∶00時(shí)(0.70 g/m3)；方案B數(shù)值最大出現(xiàn)在DOY 233 UTC 00∶00時(shí)(1.83 g/m3)，數(shù)值最小出現(xiàn)在DOY 214 UTC 12∶00時(shí)(0.54 g/m3)。方案B有近50%的試驗(yàn)RMSE小于1 g/m3，而方案A僅4組試驗(yàn)RMSE小于1 g/m3，同時(shí)方案B有90%以上的RMSE小于方案A，故單天解精度對(duì)比中方案B明顯優(yōu)于方案A。

圖5 層析時(shí)段內(nèi)兩種方案層析結(jié)果的RMSE對(duì)比Fig.5 Comparison of the RMSE of the tomography results derived from 2 projects during the tomography periods

本文以高精度和高垂直分辨率的探空水汽密度廓線信息為參考值，對(duì)整個(gè)試驗(yàn)數(shù)據(jù)集的層析結(jié)果進(jìn)行整體精度評(píng)定，表1給出了兩種方案層析結(jié)果的RMSE和MAE(mean absolute error)。由表1可知，方案B的兩種統(tǒng)計(jì)量相比于方案A均具有更低數(shù)值，其中RMSE降低了0.401 g/m3，MAE降低了0.223 g/m3，說明ANEVS的整體解算精度高于均勻分層。

表1 兩種方案的層析結(jié)果與探空數(shù)據(jù)對(duì)比的平均RMSE和平均MAETab.1 Average RMSE and average MAE of the tomography results derived from 2 projects using radiosonde data as a reference (g/m3)

圖6為試驗(yàn)時(shí)段內(nèi)，兩種方案層析水汽密度與探空數(shù)據(jù)之間的散點(diǎn)圖。如圖6所示，方案B相較于方案A具有更強(qiáng)的聚集性，同時(shí)方案B具有更高的R2、更低的RMSE，擬合線性函數(shù)的斜率也更接近1，表明方案B解算結(jié)果質(zhì)量?jī)?yōu)于方案A。對(duì)比兩方案的散點(diǎn)圖還可發(fā)現(xiàn)，方案B散點(diǎn)在20～25 g/m3區(qū)間內(nèi)的聚集程度高于方案A，這是由于均勻分層于近地區(qū)域的水汽分層粗糙，底層間距大，水汽密度跳躍殘缺，衛(wèi)星信號(hào)利用率低，大幅降低了方案A底層的解算精度，導(dǎo)致底層解算精度與中間層解算精度相當(dāng)，因此方案A的散點(diǎn)圖于20～25 g/m3區(qū)間內(nèi)的聚集程度十分接近5～20 g/m3區(qū)間，而方案B的底層水汽分層細(xì)致，信號(hào)利用率高，大幅提高了底層解算精度。下面本文將對(duì)兩種方案在不同天氣條件下和不同高度處的層析結(jié)果做進(jìn)一步分析。

圖6 兩種方案層析反演水汽密度與探空水汽密度散點(diǎn)圖Fig.6 Scatter plots of water vapor density between radiosonde and tomography derived from 2 projects

4.2 不同天氣層析結(jié)果分析

圖7為層析水汽密度廓線與探空水汽密度對(duì)比圖，其中紅色曲線和藍(lán)色曲線分別表示方案A和方案B的垂直水汽廓線信息，綠色散點(diǎn)表示探空水汽密度信息。本文選取DOY 221 UTC 00∶00時(shí)和DOY 238 UTC 00∶00時(shí)的層析水汽廓線為例與探空數(shù)據(jù)進(jìn)行驗(yàn)證分析，其中DOY 221最高溫度達(dá)34.8℃，無降水活動(dòng)發(fā)生，屬于高溫炎熱的無雨天氣；DOY 238由于受本月第二次強(qiáng)臺(tái)風(fēng)“白鹿”的影響出現(xiàn)持續(xù)大暴雨，屬于暴雨天氣。由圖7可知，方案B與探空水汽密度具有更好的符合度，尤其在3000 m以下的高度區(qū)間內(nèi)，方案B的水汽廓線能更好地反映水汽的真實(shí)情況。此外，對(duì)比兩種方案在不同天氣條件下的垂直水汽廓線可以發(fā)現(xiàn)，方案A在近地高度區(qū)間內(nèi)僅能表現(xiàn)出大氣水汽的整體變化趨勢(shì)，其原因是在進(jìn)行層析區(qū)域建模時(shí)，基于均勻分層的方案A在3000 m以下僅具有3層完整的分層區(qū)間，導(dǎo)致通過水汽層析技術(shù)反演獲得的水汽密度參數(shù)數(shù)量較少，無法表現(xiàn)出大氣水汽在惡劣天氣條件下的波動(dòng)現(xiàn)象，而基于ANEVS的方案B在3000 m以下具有8層完整的分層區(qū)間，充足的水汽密度參數(shù)使其具備了反映近地空間大氣水汽波動(dòng)現(xiàn)象的能力。

圖7 兩種方案的層析水汽廓線與探空水汽密度對(duì)比Fig.7 Comparison of tomographic water vapor profiles between 2 projects using the radiosonde as a reference

水汽密度廓線只具備表現(xiàn)天頂方向的垂直水汽密度信息的能力，并不能對(duì)三維水汽密度場(chǎng)進(jìn)行展示，如圖8所示，本節(jié)基于方案A和方案B的層析反演結(jié)果，展示了DOY 238 UTC 00∶00時(shí)三維水汽密度場(chǎng)的可視化結(jié)果。方案A與方案B相比，方案A中下邊界處層與層之間的色階變化明顯，上邊界處的色階變化則較小，表明均勻分層反演的三維水汽密度信息具有下邊界處匱乏而上邊界處冗余的缺點(diǎn)，對(duì)比圖中方案B反演的三維水汽場(chǎng)可發(fā)現(xiàn)，ANEVS彌補(bǔ)了均勻分層的不足，其色階變化豐富，具有更精確的水汽密度信息反演能力。

圖8 兩種方案的三維水汽密度場(chǎng)可視化Fig.8 Visualization of three dimensional water vapor density field derived from 2 projects

表2為DOY 221和DOY 238不同方案層析水汽密度與探空水汽密度的相關(guān)系數(shù)。在無雨日(DOY 221)和降雨日(DOY 238)的對(duì)比中，兩種方案均在無雨天氣情況下具有更高的相關(guān)系數(shù)，說明惡劣天氣一定程度上影響了層析結(jié)果的質(zhì)量；在兩種方案之間的對(duì)比中，不同天氣條件下方案B的相關(guān)系數(shù)均高于方案A，表明ANEVS相比于傳統(tǒng)均勻分層具有更優(yōu)質(zhì)的大氣水汽反演效果；對(duì)比兩種天氣條件下不同方案間的相關(guān)系數(shù)之差可發(fā)現(xiàn)，在水汽不穩(wěn)定的惡劣天氣，方案B對(duì)于相關(guān)系數(shù)的提高量大于水汽穩(wěn)定的無雨天氣，說明方案B在惡劣天氣條件下能更好地優(yōu)化大氣水汽的反演效果。

4.3 不同高度層析結(jié)果分析

圖9給出了兩種方案的層析結(jié)果在每1層上與探空數(shù)據(jù)對(duì)比的RMSE。從圖中可看出，兩種分層方案的各層RMSE均呈現(xiàn)先增大再減小的變化規(guī)律，在方案B中所有分層的RMSE均小于1.5 g/m3，而方案A中1～5層的RMSE大于1.5 g/m3，其中第2層的RMSE大于2.5 g/m3，表明方案A于低海拔處的層析解算結(jié)果質(zhì)量較低。從兩種方案的各層RMSE對(duì)比中還可發(fā)現(xiàn)，方案A中9～11層的RMSE小于1 g/m3，而方案B中對(duì)應(yīng)層數(shù)的RMSE大于1 g/m3，產(chǎn)生此現(xiàn)象的原因是兩種方案在分層數(shù)相同的情況下，每1層所對(duì)應(yīng)的實(shí)際高度區(qū)間不同，分層層數(shù)越大，兩個(gè)方案之間的實(shí)際高度區(qū)間相差越大。由4.2節(jié)可知，分層數(shù)9～11層于方案A中實(shí)際對(duì)應(yīng)6628～9113.5 m，于方案B中其實(shí)際對(duì)應(yīng)2688～4982 m，而6628～9113.5 m區(qū)間內(nèi)的實(shí)際水汽含量相較于2688～4982 m區(qū)間內(nèi)的小很多，導(dǎo)致方案A的RMSE在9～11層小于方案B。因此，為了進(jìn)一步對(duì)比兩種方案的垂直精度，圖9還給出了方案A和方案B在各高度處的RMSE對(duì)比圖。由圖9可知，方案A于4000 m以下的解算精度明顯低于方案B，這是由于方案A低海拔處分層粗糙，不符合大氣水汽的垂直分布特征，衛(wèi)星信號(hào)利用率低下，大幅降低了層析解算精度；4000 m以上兩種方案精度相當(dāng)，即使方案A于高海拔處分層更精細(xì)，但由于衛(wèi)星信號(hào)分布稀疏、水汽稀薄，解算精度并未得到實(shí)質(zhì)性提高。方案B使用ANEVS由于考慮到水汽的垂直變化特征，以水汽密度區(qū)間為分層準(zhǔn)則，該方法獲得的非均勻高度區(qū)間使各層水汽密度變化控制在一個(gè)較小的范圍內(nèi)，同時(shí)這種非均勻高度區(qū)間還適應(yīng)了衛(wèi)星信號(hào)呈倒錐形發(fā)散的特征，在保證了高海拔處解算絕對(duì)誤差較小的同時(shí)，大幅提高了低海拔處的解算精度。

圖9 兩種方案層析結(jié)果與探空數(shù)據(jù)對(duì)比的各分層RMSEFig.9 RMSE at each layer of the tomographic results derived from 2 projects using the radiosonde as a reference

大氣水汽整體表現(xiàn)出高層稀薄、底層充沛的特點(diǎn)，這種現(xiàn)象使高層水汽密度的絕對(duì)誤差相較于其他位置小很多，再加上傳統(tǒng)均勻分層將過多的層次劃分至絕對(duì)誤差本就較小的高層區(qū)域內(nèi)，因此使用與絕對(duì)誤差相關(guān)的統(tǒng)計(jì)量評(píng)定解算精度可能造成誤判，對(duì)此本文額外使用相對(duì)誤差(relative error，RE)對(duì)層析結(jié)果進(jìn)行評(píng)定。相對(duì)誤差表示為絕對(duì)誤差與真值的比值，計(jì)算公式為

(15)

式中，RE表示相對(duì)誤差(%)；DT表示水汽層析反演水汽密度值，單位為g/m3；DRS表示利用探空數(shù)據(jù)計(jì)算得到的水汽密度值，單位為g/m3。

相對(duì)誤差能更好地考慮真值本身的數(shù)量級(jí)，用百分比的形式對(duì)誤差進(jìn)行歸一化處理，因此在研究不同垂直分層方法之間解算精度時(shí)使用相對(duì)誤差具有更好的真實(shí)性和可比性。如圖10所示，兩種方案的各層平均相對(duì)誤差(mean relative error，MRE)均隨高度增加而增加，表明解算精度隨高度增加而遞減。在2500 m以下的高度區(qū)間內(nèi)，兩種方案的MRE均小于20%，解算精度較高，ANEVS的MRE在2000 m以下小于10%，在1000 m以下小于5%，精度可達(dá)傳統(tǒng)均勻分層的2倍以上，具有更高的解算精度；在2500～4500 m的高度區(qū)間內(nèi)，兩種方案的解算精度同時(shí)隨高度增加而降低，但方案B的精度仍稍優(yōu)于方案A。由圖10還可以發(fā)現(xiàn)，方案A只有1～7層MRE小于100%，8～13層MRE均大于100%，說明傳統(tǒng)均勻分層并未切實(shí)提高高層區(qū)域的解算精度；方案B僅最頂層MRE超過100%，倒數(shù)第2層MRE接近40%，其余層MRE數(shù)值更小，ANEVS細(xì)致的低層層次劃分適應(yīng)了大氣水汽的分布特征，充分高效地利用了低海拔密集的衛(wèi)星信號(hào)觀測(cè)值，水汽層析解算效率遠(yuǎn)高于方案A。綜上，在以MRE為依據(jù)的解算精度對(duì)比上，ANEVS相比于傳統(tǒng)均勻分層有很大提高，該方法能將低海拔水汽進(jìn)行細(xì)分，控制各層水汽密度差異，提高衛(wèi)星信號(hào)的利用率，將高海拔處密度相近的水汽歸并為同一層次，減少了高層水汽解算的多余工作，在自適應(yīng)地進(jìn)行大氣水汽垂直分層的同時(shí)，克服了傳統(tǒng)均勻分層方案的缺點(diǎn)。

圖10 兩種方案層析結(jié)果與探空數(shù)據(jù)對(duì)比的各層MREFig.10 MRE at each layer of the tomographic results derived from 2 projects using the radiosonde as a reference

5 結(jié) 語

本文提出了一種基于大氣水汽垂直分布特征的自適應(yīng)非均勻指數(shù)垂直分層方法(ANEVS)。ANEVS以水汽密度區(qū)間和水汽垂直分布特征函數(shù)反推垂直分層高度區(qū)間，再通過迭代自適應(yīng)地調(diào)整高度區(qū)間，以縮小各層之間的水汽密度差異，使分層高度區(qū)間滿足實(shí)際水汽變化規(guī)律，實(shí)現(xiàn)對(duì)任意給定層析區(qū)域進(jìn)行任意分層層數(shù)的自適應(yīng)非均勻分層，該非均勻分層方法在具備科學(xué)可靠分層依據(jù)的同時(shí)又具有方便快捷且自適應(yīng)的優(yōu)勢(shì)。本文利用香港CORS網(wǎng)2019年8月的GNSS觀測(cè)數(shù)據(jù)對(duì)比分析了ANEVS與傳統(tǒng)均勻分層的層析結(jié)果的質(zhì)量及精度，得到以下結(jié)論。

(1)ANEVS在層析結(jié)果的精度與質(zhì)量上均優(yōu)于傳統(tǒng)均勻分層。在散點(diǎn)圖的對(duì)比上，ANEVS表現(xiàn)出更強(qiáng)的聚集性，同時(shí)ANEVS散點(diǎn)圖具有更高的確定系數(shù)和更低的均方根誤差，ANEVS層析結(jié)果的均方根誤差和平均絕對(duì)誤差相較于傳統(tǒng)均勻分層分別降低0.401 g/m3和0.223 g/m3。

(2)在不同天氣條件下，ANEVS水汽廓線能更好地反映水汽的垂直擾動(dòng)變化。不同天氣下ANEVS水汽密度與探空水汽密度的相關(guān)系數(shù)均高于傳統(tǒng)均勻分層，在暴雨天氣下兩種方法的相關(guān)系數(shù)之差大于無雨天氣，表明ANEVS在惡劣天氣條件下對(duì)大氣水汽反演效果的優(yōu)化更加顯著。

(3)在不同高度的層析結(jié)果對(duì)比中，ANEVS對(duì)低層高度區(qū)間內(nèi)層析結(jié)果精度的提高更加顯著。在分層數(shù)同為13層的條件下，傳統(tǒng)均勻分層的平均相對(duì)誤差有6層高于100%，僅1～7層小于50%，ANEVS的平均相對(duì)誤差僅最頂層1層高于100%，1～9層小于20%，ANEVS的均方根誤差于4000 m以下的高度區(qū)間內(nèi)顯著減小，表明ANEVS大幅提高了各層層析結(jié)果的解算質(zhì)量，其中低層層析結(jié)果質(zhì)量的提升尤為顯著。

筆者后續(xù)的研究與工作內(nèi)容主要包括兩方面：一方面是調(diào)整并擴(kuò)大層析試驗(yàn)區(qū)域，以便引入復(fù)數(shù)個(gè)探空氣象站的數(shù)據(jù)，對(duì)大氣水汽垂直分布特征模型的有效覆蓋范圍進(jìn)行定量分析；另一方面，基于對(duì)先驗(yàn)大氣水汽密度進(jìn)行長(zhǎng)時(shí)間跨度的統(tǒng)計(jì)學(xué)分析，研究并總結(jié)出一種能夠適用于不同月份的表示大氣水汽垂直分布特征的建模準(zhǔn)則。

致謝：特別感謝香港地政總署測(cè)繪處提供的香港衛(wèi)星定位參考站網(wǎng)絡(luò)的GNSS觀測(cè)數(shù)據(jù)(https:∥www.geodetic.gov.hk/en/satref/satref.htm)，香港King 's Park天文臺(tái)提供的京士柏氣象站的無線電探空數(shù)據(jù)(http:∥weather.uwyo.edu/upperair/sounding.html)，麻省理工學(xué)院地球大氣和行星科學(xué)系提供的GAMIT/GLOBK軟件。