顧及邊界信號及垂直約束的GNSS水汽層析方法

何秀鳳，詹 偉，施宏凱

河海大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 211100

精確的大氣水汽三維時(shí)空分布信息對于提升中尺度數(shù)值氣象預(yù)報(bào)精度具有非常重要的意義[1]。當(dāng)前對大氣水汽進(jìn)行探測的常規(guī)技術(shù)(無線電探空、氣象衛(wèi)星等)存在低時(shí)空分辨率、高成本、觀測精度受天氣影響大、觀測范圍局限等問題[2-3]。全球?qū)Ш叫l(wèi)星系統(tǒng)(GNSS)最初被應(yīng)用于區(qū)域大氣可降水量[4-6]的估計(jì)，此后為了進(jìn)一步得到區(qū)域內(nèi)水汽的三維分布，GNSS三維水汽層析技術(shù)[7]興起，憑借其低成本、高時(shí)空分辨率、全天候的特點(diǎn)被廣泛應(yīng)用。但傳統(tǒng)算法[8-9]只考慮從層析區(qū)域頂端射出的信號，忽略了側(cè)面穿刺的信號，層析方法存在觀測數(shù)據(jù)利用率低、觀測分布不均等問題。有效利用高度角較低的邊界信號提供的信息可以增加層析區(qū)域內(nèi)網(wǎng)格尤其是底層及邊緣網(wǎng)格被穿過的概率、減小層析網(wǎng)格空格率[10]、提高觀測數(shù)據(jù)利用率。

為解決邊界信號未能有效利用的問題，文獻(xiàn)[11—13]建立了一種基于水汽密度比例因子的三維水汽層析算法；文獻(xiàn)[14]通過縮小層析水平范圍構(gòu)建層頂入射信號上的水汽含量與范圍縮小后的層析邊界高程間的比例關(guān)系，用于恢復(fù)邊界信號的有效水汽信息；文獻(xiàn)[15—16]附加一個(gè)輔助層析區(qū)域，使從研究區(qū)域側(cè)面入射的信號也可被有效利用。以上3種模型在一定程度上克服了觀測數(shù)據(jù)利用率較低的局限性，取得了較高的精度，但是均在單一導(dǎo)航系統(tǒng)中得以實(shí)現(xiàn)，在多模組合系統(tǒng)水汽層析中尚未有所研究，觀測數(shù)據(jù)利用率有待進(jìn)一步提高。在此背景下，文獻(xiàn)[10]對以上研究進(jìn)行了詳細(xì)分析，通過層析時(shí)間前3 d的探空信息擬合出不同高程對應(yīng)的大氣可降水量(precipitable water vapor,PWV)與高程的函數(shù)關(guān)系，從PWV的角度恢復(fù)邊界信號水汽信息，在多模組合系統(tǒng)中進(jìn)行水汽層析試驗(yàn)，并證實(shí)了引入多系統(tǒng)信號及引入邊界信號均能獲取更好的水汽分布結(jié)果。

文獻(xiàn)[10]基于擬合函數(shù)建立的邊界信號的精度相對于層頂信號存在較大差別，有必要削弱這部分誤差對層析結(jié)果的影響。與此同時(shí)，以上研究中的層析模型均是根據(jù)垂直方向上水汽分布呈指數(shù)遞減規(guī)律[17]建立上下層的關(guān)系，得到的傳統(tǒng)垂直約束方程與層析區(qū)域?qū)嶋H水汽分布符合程度較低，采用傳統(tǒng)垂直約束解算的結(jié)果精度較差，在底層區(qū)域表現(xiàn)得較為明顯。為了進(jìn)一步提高三維水汽層析的精度，有必要對上述兩個(gè)問題加以考慮。因此，本文基于短期探空信息擬合函數(shù)建立了層析區(qū)域內(nèi)水汽分布的區(qū)域性垂直約束，采用了一種顧及邊界信號的GNSS水汽層析方法，設(shè)計(jì)和實(shí)現(xiàn)了融合邊界信號的多模水汽層析試驗(yàn)，采用無線電探空數(shù)據(jù)進(jìn)行精度驗(yàn)證，詳細(xì)對比分析了邊界信號及本文建立的垂直約束對水汽層析結(jié)果的改善程度。

1 基于探空信息擬合函數(shù)的三維水汽層析算法

1.1 基本觀測方程

傾斜路徑水汽含量(slant water vapor,SWV)指的是GNSS信號傳播方向上的水汽柱內(nèi)所含水汽轉(zhuǎn)化為單位面積上的液態(tài)水高度，單位為mm，其表達(dá)式為

(1)

式中，S為SWV；L為衛(wèi)星到接收機(jī)的信號路徑；ρwater表示水汽密度，單位為g/m3；dL表示單位截面積上的體長度，單位為km。將式(1)方程離散化為線性形式的觀測方程

S=∑ijk(dijk·ρijk)=A·X

(2)

式中，dijk為信號ρ在網(wǎng)格(i、j、k)內(nèi)的截距，單位為km；ρijk為(i、j、k)網(wǎng)格內(nèi)的水汽密度；A表示斜路徑信號在層析網(wǎng)格內(nèi)穿過的截距向量，未穿過的網(wǎng)格截距為零；X為層析網(wǎng)格內(nèi)部各個(gè)網(wǎng)格的水汽密度值。根據(jù)式(2)將所有GNSS觀測表達(dá)式表示為矩陣形式

Am×n·Xn×1=Sm×1

(3)

式中，Sm×1=[S1S2…Sm]T表示所有GNSS斜路徑信號上觀測值向量；Am×n=[A1A2…Am]T為對應(yīng)信號在網(wǎng)格內(nèi)的截距矩陣；Xn×1為各個(gè)網(wǎng)格的水汽密度未知量；m表示衛(wèi)星信號的條數(shù)；n為網(wǎng)格個(gè)數(shù)。

1.2 邊界信號觀測方程

針對傳統(tǒng)模型中未能有效利用許多側(cè)邊穿刺信號的問題，文獻(xiàn)[10]提出直接從PWV與高程的函數(shù)關(guān)系恢復(fù)邊界信號水汽值的方法。其思路是：PWV作為求解SWV的關(guān)鍵元素，與高程存在函數(shù)關(guān)系，而與梯度改正項(xiàng)無關(guān)，因而利用五階多項(xiàng)式擬合PWV與高程的相互關(guān)系，直接恢復(fù)邊界信號的SWV。本文考慮到水汽在垂直方向上呈指數(shù)遞減的規(guī)律，利用指數(shù)函數(shù)擬合兩者的關(guān)系以恢復(fù)邊界信號信息，并為層析高度的選取提供參考。其中，PWV與高程相互關(guān)系擬合過程如下。

(1)利用層析區(qū)域內(nèi)探空站前3 d的探空數(shù)據(jù)，在垂直方向上計(jì)算出采集點(diǎn)累計(jì)PWV占整個(gè)垂直方向上PWV的比例因子γ。

(2)由步驟(1)計(jì)算的γ及其對應(yīng)的高程H進(jìn)行指數(shù)函數(shù)擬合，建立函數(shù)關(guān)系式γ=a·eb·H+c·ed·H。

(3)通過步驟(2)得到的擬合函數(shù)，只需要代入每條邊界信號穿刺點(diǎn)處的高程H，便可由測站處PWV計(jì)算得到對應(yīng)的位于層析區(qū)域內(nèi)的有效PWVH=PWV·γ。然后采用SWV計(jì)算公式[18]計(jì)算出邊界信號水汽值SWVH

SWVH=mw(θ)·PWVH+Π·{mΔ(θ)·

(GN·cosφ+GE·sinφ)+ε}

(4)

式中，mw(θ)、mΔ(θ)分別為濕映射函數(shù)和梯度映射函數(shù)，與衛(wèi)星高度角θ有關(guān)；Π為水汽轉(zhuǎn)換系數(shù)；GN、GE分別為南北方向和東西方向梯度參數(shù)；φ為方位角。綜上建立起邊界入射信號的觀測方程，矩陣形式如下

AHm×n·Xn×1=SHm×1

(5)

式中，SHm×1為邊界信號位于層析區(qū)域內(nèi)的有效水汽值；AHm×n為邊界信號的系數(shù)矩陣；m為信號條數(shù)。

1.3 水平約束與垂直約束

GNSS觀測方程系數(shù)矩陣是一個(gè)絕大多數(shù)為零的稀疏矩陣，通過引入額外的水平和垂直約束條件[19]對觀測方程進(jìn)行改善。

水平約束條件根據(jù)同層網(wǎng)格之間的相關(guān)性，基于高斯加權(quán)法[11]對每一個(gè)網(wǎng)格建立約束條件，構(gòu)成如下水平約束方程

H·X=0

(6)

式中，H為水平約束系數(shù)矩陣。

文獻(xiàn)[20]通過大量試驗(yàn)指出，對流層低層的垂直分辨率相對較高，采用垂直不均勻分層方法得到的大氣水汽更加精確，垂直不均勻分層的最佳高度為500～1200 m。傳統(tǒng)方法中垂直約束根據(jù)垂直方向上水汽分布呈指數(shù)遞減規(guī)律建立上下層的關(guān)系

xk+1=xk·e(-Z/H)

(7)

式中，xk+1、xk分別為上下兩層的水汽密度參數(shù)；Z為上下層的高程差；H為水汽標(biāo)高，通常取值1～2 km。

本文垂直約束由1.2小節(jié)中得到的擬合函數(shù)構(gòu)建上下兩層之間的函數(shù)關(guān)系

xk+1=α·xk

(8)

V·X=0

(9)

式中，V為垂直約束方程系數(shù)矩陣。

1.4 層析方程組解算

綜合上述內(nèi)容，引入GNSS多模觀測數(shù)據(jù)得到顧及邊界信號的層析方程組

(10)

式中，m=g+r+e+c；Hm=Hg+Hr+He+Hc；m、g、r、e、c分別代表GNSS組合和G、R、E、C各單系統(tǒng)的層頂入射信號數(shù)目；Hm、Hg、Hr、He、Hc分別代表GNSS組合和G、R、E、C各單系統(tǒng)的邊界入射信號數(shù)目，即

添加約束條件后的層析方程組仍然是一個(gè)絕大部分元素為零的稀疏矩陣，如何對方程組進(jìn)行解算將影響層析結(jié)果的準(zhǔn)確性及穩(wěn)定性，一般采用代數(shù)重構(gòu)法[21-25](algebraic reconstruction techniques,ARTs)求解水汽參數(shù)，可以避免系數(shù)矩陣求逆的運(yùn)算。本文采用加型ART進(jìn)行方程組解算，其公式如下

(11)

式中，k表示第k次迭代；xk為第k次迭代的結(jié)果；Ai表示觀測矩陣的第i行；Si為對應(yīng)的SWVi；〈·,·〉表示向量內(nèi)積；λ為松弛因子。

2 GNSS三維水汽層析

2.1 試驗(yàn)區(qū)域及數(shù)據(jù)

選擇香港地區(qū)作為試驗(yàn)區(qū)域，獲取了香港CORS的10個(gè)站點(diǎn)2016年DOY 129—DOY 135一周的GNSS觀測數(shù)據(jù)和氣象數(shù)據(jù)，采用CODE分析中心提供的精密星歷、鐘差、軌道等產(chǎn)品，基于數(shù)據(jù)處理軟件RTKLIB進(jìn)行GPS/GLONASS/Galileo/BDS四系統(tǒng)組合精密單點(diǎn)定位(precise point positioning,PPP)解算，估算對流層延遲參數(shù)，數(shù)據(jù)處理策略見表1。

表1 PPP數(shù)據(jù)處理策略Tab.1 PPP calculation strategy

研究區(qū)域網(wǎng)格劃分：緯度范圍為22.16°N—22.52°N，經(jīng)度范圍為113.90°E—114.35°E，經(jīng)緯度網(wǎng)格密度為8×8均勻劃分；垂直方向上由探空信息擬合函數(shù)計(jì)算得到8000 m高度(DOY 129—DOY 135)PWV比例因子均值為0.982，最小值為0.973，8000 m高度以上分布的水汽均值低于0.9 mm，最大低于1.6 mm，因而將層析高度定為8000 m。采用不均勻分層(0、500、1000、1700、2400、3100、4000、4900、5800、6900、8000，單位為m)，總計(jì)640個(gè)網(wǎng)格。利用KingsPark探空站無線電探空數(shù)據(jù)作為水汽層析結(jié)果的驗(yàn)證數(shù)據(jù)。研究區(qū)域及測站分布如圖1所示。

圖1 層析區(qū)域及GNSS測站分布Fig.1 Tomography area and distribution of GNSS stations

2.2 探空信息擬合函數(shù)精度分析

圖2為2016年DOY 129—DOY 134根據(jù)探空站數(shù)據(jù)計(jì)算的PWV比例因子和高程相互關(guān)系的擬合結(jié)果。由圖2可知：DOY 129—DOY 131 3 d采樣點(diǎn)數(shù)據(jù)分布相對密集，對應(yīng)的擬合結(jié)果均方根誤差(RMSE)分別為0.019 2、0.017 3、0.021 1；DOY 132—DOY 134 3 d則相對分散一些，RMSE分別為0.039 1、0.049 9、0.052 8。圖2結(jié)果表明，基于探空信息得到的擬合函數(shù)存在一定的模型誤差，表現(xiàn)在PWVH上1～3 mm，由此模型計(jì)算的邊界信號的SWV觀測值相對于層頂信號的SWV精度要低一些，因此，有必要在層析方程組解算時(shí)采取合理的方法對這部分誤差進(jìn)行削弱。

2.3 層析區(qū)域內(nèi)信號數(shù)量與空格率分析

為研究邊界信號對層析模型的改善程度，本文首先統(tǒng)計(jì)了僅考慮層析區(qū)域天頂信號(圖3(a)，以GNSS表示)和顧及邊界信號(圖3(b)，以GNSS+表示)的兩種層析模型在30 min內(nèi)的有效GNSS信號數(shù)量和空格率，其中數(shù)據(jù)采樣間隔為1 min。

圖2 DOY 129—DOY 134指數(shù)函數(shù)擬合曲線Fig.2 Curve fitted by exponential function of DOY 129—DOY 134

圖3 兩種層析模型Fig.3 Diagrammatic drawing of two tomographic models

融合邊界信號后層析區(qū)域內(nèi)的觀測量增多，同時(shí)低高度角的邊界信號使得層析區(qū)域底層及邊緣網(wǎng)格被穿過的機(jī)率大大提高。為此，本文對DOY 129—DOY 135每天UTC 12 h層析模型的觀測量、網(wǎng)格空格率及邊緣網(wǎng)格(層析區(qū)域側(cè)面最外圍網(wǎng)格，不包括頂層)空格率進(jìn)行了對比分析。其中層析區(qū)域內(nèi)網(wǎng)格總量為640、邊緣網(wǎng)格為198，統(tǒng)計(jì)結(jié)果見圖4和表2。

圖4和表2的結(jié)果表明：顧及邊界信號的GNSS+模型在觀測量、網(wǎng)格空格率，以及邊緣網(wǎng)格空格率3個(gè)方面上相對于舍棄邊界信號的GNSS模型都存在優(yōu)勢。在層析解算時(shí)刻：GNSS方案能夠利用的平均有效觀測量只有4637，而GNSS+為7045，觀測量提高了51.9%，并且GNSS+模型的最低觀測量都要高于GNSS模型的最高觀測量；GNSS+方案的平均空格率為14.3%，相對于GNSS的26.1%降低了11.8%；邊緣空格率則表現(xiàn)得更為明顯，相對于GNSS方案47.1%的平均空格率，GNSS+降低了24.1%，只有23.0%。

2.4 三維層析結(jié)果精度分析

本文利用2016年DOY 129—DOY 135一周的數(shù)據(jù)進(jìn)行三維水汽層析精度分析，Kingspark探空站位于層析區(qū)域內(nèi)，并且其配備的無線電探空儀能夠提供垂直方向上精確的水汽密度廓線圖，因而將該探空站所在網(wǎng)格對應(yīng)的層析垂直廓線與探空數(shù)據(jù)進(jìn)行對比。本文采用平均絕對偏差(MAE)和均方根誤差(RMS)作為評定層析方法的指標(biāo)，共設(shè)計(jì)以下4種方案進(jìn)行試驗(yàn)。

圖4 有效信號數(shù)、空格率、邊緣空格率統(tǒng)計(jì)Fig.4 Statistical results of effective signal number,space rate and edge space rate

表2 有效信號數(shù)、空格率、邊緣空格率統(tǒng)計(jì)Tab.2 Statistical results of effective signal number,space rate and edge space rate

方案1：僅考慮層析區(qū)域天頂信號的層析模型(GNSS)。

方案2：顧及邊界信號的層析模型(GNSS+)。

方案3：顧及邊界信號但不考慮PWV擬合函數(shù)模型誤差的層析模型(GNSS-)。

方案4：顧及邊界信號但使用傳統(tǒng)垂直約束的層析模型(TRA)。

為了對比分析4種方案的層析效果，本文首先計(jì)算了DOY 129—DOY 134每天UTC 12 h各種方案解算的無線電探空站所在位置不同高度上的水汽密度，并分別與無線電探空儀計(jì)算的水汽密度對比。各方案與無線電探空儀對比的垂直廓線結(jié)果如圖5所示。

圖5表明，各方案在層析時(shí)刻的水汽結(jié)果具有相同的變化趨勢，與探空水汽結(jié)果的吻合程度均較高，但是各方案的反演結(jié)果在底層區(qū)域有所差異，使用傳統(tǒng)垂直約束的TRA模型反演的底層水汽密度值易出現(xiàn)偏大的情況，這是由于層析區(qū)域底層水汽密度并非指數(shù)分布，采用按照指數(shù)遞減規(guī)律的傳統(tǒng)垂直約束方程在底層會導(dǎo)致較大偏差。為了對層析結(jié)果和探空數(shù)據(jù)計(jì)算結(jié)果進(jìn)行直接對比，本文統(tǒng)計(jì)了DOY 129—DOY 135每天UTC 12 h各種方案的層析結(jié)果，利用層析結(jié)果計(jì)算出無線電探空站所在位置的水汽密度，分別與無線電探空儀數(shù)據(jù)計(jì)算的結(jié)果進(jìn)行對比。圖6為各方案統(tǒng)計(jì)的每天的MAE和RMS，表3為對應(yīng)的統(tǒng)計(jì)結(jié)果。

圖5 DOY 129—DOY 134(UTC 12)時(shí)的垂直水汽廓線分布Fig.5 Vertical water vapor profiles distribution at UTC 12 h in DOY 129—DOY 134

圖6 4種方案(DOY 129—DOY 135)層析結(jié)果與無線電探空儀對比誤差統(tǒng)計(jì)Fig.6 The statistical result of difference between different tomographic results and radiosonde for DOY 129—DOY 135

表3 4種方案(DOY 129—DOY 135)層析結(jié)果與無線電探空儀對比誤差統(tǒng)計(jì)Tab.3 The statistical result of difference between different tomographic results and radiosonde for DOY 129—DOY 135 g/m3

圖6表明，GNSS方案和TRA方案整體上層析結(jié)果較差，MAE和RMS偏大，GNSS+方案結(jié)果最優(yōu)，GNSS-方案略低于GNSS+方案。由表4數(shù)據(jù)可知，GNSS+方案的層析結(jié)果平均MAE和RMS最優(yōu)，分別為1.45 g/m3和1.82 g/m3，相對于GNSS-方案略有降低，分別為1.4%和1.6%；相對于GNSS方案和TRA方案有較大降低，前者分別降低9.4%和12.1%，后者分別降低5.2%和5.7%。這說明本文采用的顧及邊界信號并考慮PWV擬合函數(shù)模型誤差的GNSS+層析模型精度最高，相對于不考慮PWV擬合函數(shù)模型誤差的GNSS-模型精度略有提升，相對于只考慮層頂信號的GNSS模型和使用傳統(tǒng)垂直約束的TRA模型這兩個(gè)傳統(tǒng)的層析模型，精度存在較大提高。

為了比較不同方案下層析結(jié)果的精確性和可靠性，本文將DOY 129—DOY 135一周數(shù)據(jù)下4種方案解算的探空站上廓線水汽密度與無線電探空數(shù)據(jù)進(jìn)行相關(guān)性分析，結(jié)果如圖7所示。

圖7 4種方案下層析結(jié)果與探空結(jié)果的線性相關(guān)程度Fig.7 Linear correlation between the results of tomography and the results of sounding information of four schemes

圖7表明，GNSS+層析模型的相關(guān)系數(shù)最高，達(dá)到0.946。GNSS-方案相關(guān)系數(shù)略低于GNSS+。GNSS和TRA兩種傳統(tǒng)模型相關(guān)系數(shù)較低，GNSS最低。相關(guān)系數(shù)分析結(jié)果與上述MAE和RMS分析結(jié)果相一致，GNSS+方案結(jié)果最優(yōu)，精度最高。此外，TRA模型在底層區(qū)域解算的水汽密度值(圖中方框內(nèi)區(qū)域)誤差較大，這與TRA模型在底層區(qū)域水汽密度解算結(jié)果易偏大的特點(diǎn)相一致。

3 結(jié) 論

本文以香港CORS網(wǎng)2016年5月8日至2016年5月14日的觀測數(shù)據(jù)和無線電探空數(shù)據(jù)，對融合邊界信號、基于探空信息擬合函數(shù)建立垂直約束的多模層析方法進(jìn)行了試驗(yàn)論證，詳細(xì)分析了邊界信號及垂直約束方法對水汽反演結(jié)果的優(yōu)化效果。研究結(jié)果表明：邊界信號的加入使得觀測量提高了51.9%，空格率尤其是邊緣和底層空格率大為下降，層析結(jié)果的平均MAE和RMS分別降低了9.4%、12.1%；通過搜索松弛因子削弱了PWVH擬合函數(shù)誤差對層析結(jié)果的影響；垂直約束方法的改進(jìn)改善了水汽反演結(jié)果，能夠改善基于指數(shù)遞減特性的傳統(tǒng)垂直約束在底層區(qū)域反演結(jié)果易偏大的不足，平均MAE和RMS分別降低了5.2%、5.7%。本文采用的GNSS+方案的層析結(jié)果相較其他方案能夠獲得更高精度的三維水汽分布信息。