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

何秀鳳，詹 偉，施宏凱

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

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

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

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

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

1.1 基本觀測方程

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

(1)

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

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

(2)

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

Am×n·Xn×1=Sm×1

(3)

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

1.2 邊界信號觀測方程

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

(1)利用層析區域內探空站前3 d的探空數據，在垂直方向上計算出采集點累計PWV占整個垂直方向上PWV的比例因子γ。

(2)由步驟(1)計算的γ及其對應的高程H進行指數函數擬合，建立函數關系式γ=a·eb·H+c·ed·H。

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

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

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

(4)

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

AHm×n·Xn×1=SHm×1

(5)

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

1.3 水平約束與垂直約束

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

水平約束條件根據同層網格之間的相關性，基于高斯加權法[11]對每一個網格建立約束條件，構成如下水平約束方程

H·X=0

(6)

式中，H為水平約束系數矩陣。

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

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

(7)

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

本文垂直約束由1.2小節中得到的擬合函數構建上下兩層之間的函數關系

xk+1=α·xk

(8)

V·X=0

(9)

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

1.4 層析方程組解算

綜合上述內容，引入GNSS多模觀測數據得到顧及邊界信號的層析方程組

(10)

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

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

(11)

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

2 GNSS三維水汽層析

2.1 試驗區域及數據

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

表1 PPP數據處理策略Tab.1 PPP calculation strategy

研究區域網格劃分：緯度范圍為22.16°N—22.52°N，經度范圍為113.90°E—114.35°E，經緯度網格密度為8×8均勻劃分；垂直方向上由探空信息擬合函數計算得到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)，總計640個網格。利用KingsPark探空站無線電探空數據作為水汽層析結果的驗證數據。研究區域及測站分布如圖1所示。

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

2.2 探空信息擬合函數精度分析

圖2為2016年DOY 129—DOY 134根據探空站數據計算的PWV比例因子和高程相互關系的擬合結果。由圖2可知：DOY 129—DOY 131 3 d采樣點數據分布相對密集，對應的擬合結果均方根誤差(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結果表明，基于探空信息得到的擬合函數存在一定的模型誤差，表現在PWVH上1～3 mm，由此模型計算的邊界信號的SWV觀測值相對于層頂信號的SWV精度要低一些，因此，有必要在層析方程組解算時采取合理的方法對這部分誤差進行削弱。

2.3 層析區域內信號數量與空格率分析

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

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

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

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

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

2.4 三維層析結果精度分析

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

圖4 有效信號數、空格率、邊緣空格率統計Fig.4 Statistical results of effective signal number,space rate and edge space rate

表2 有效信號數、空格率、邊緣空格率統計Tab.2 Statistical results of effective signal number,space rate and edge space rate

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

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

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

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

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

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

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

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

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

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

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

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

圖7表明，GNSS+層析模型的相關系數最高，達到0.946。GNSS-方案相關系數略低于GNSS+。GNSS和TRA兩種傳統模型相關系數較低，GNSS最低。相關系數分析結果與上述MAE和RMS分析結果相一致，GNSS+方案結果最優，精度最高。此外，TRA模型在底層區域解算的水汽密度值(圖中方框內區域)誤差較大，這與TRA模型在底層區域水汽密度解算結果易偏大的特點相一致。

3 結 論

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