趙慶志，姚宜斌，姚頑強，陳 鵬，吳滿意

1. 西安科技大學(xué)測繪科學(xué)與技術(shù)學(xué)院，陜西 西安 710054； 2. 武漢大學(xué)測繪學(xué)院，湖北 武漢 430079； 3. 國家測繪地理信息局第一地形測量隊，陜西 西安 710054

地基GNSS水汽層析技術(shù)是獲取中小尺度三維水汽時空分布的重要方法之一，對于水汽的變化，特別是垂向水汽的流動具有重要的指示作用[1-2]。文獻[3]首次提出了利用區(qū)域觀測網(wǎng)中20個GPS觀測站數(shù)據(jù)重構(gòu)對流層水汽結(jié)構(gòu)的概念。文獻[4]實現(xiàn)了基于GPS層析技術(shù)的區(qū)域三維水汽信息獲取。

在國內(nèi)，文獻[5]率先對三維水汽層析問題進行了詳細(xì)、系統(tǒng)的介紹，并聯(lián)合數(shù)值模式的預(yù)報濕度場作為背景場對上海區(qū)域的三維水汽進行反演。文獻[6]提出了將高斯約束作為Kalman濾波解算初始狀態(tài)變量，用于反演濕折射率廓線的方法。文獻[1]對水汽層析中各種約束條件進行分析，發(fā)現(xiàn)水平約束方程的權(quán)陣和測站高差在不同情況下對層析結(jié)果影響不同。文獻[7]發(fā)展了基于Kalman濾波的水汽層析算法，該方法能夠高效、穩(wěn)定地層析出大氣水汽的垂直結(jié)構(gòu)，使得層析結(jié)果更加忠實于實際的水汽分布。文獻[8]通過結(jié)合數(shù)值積分參數(shù)化方法，發(fā)展了一種顧及地球曲率的三維水汽層析方法。文獻[2]對代數(shù)重構(gòu)算法在水汽層析中的應(yīng)用問題進行了詳細(xì)討論，并給出了最優(yōu)松弛因子的黃金分割搜索法和確定終止條件的NCP規(guī)則。文獻[9]提出了基于代數(shù)重構(gòu)算法的GNSS水汽層析方法，該方法能夠節(jié)省計算機內(nèi)存且穩(wěn)定度高。文獻[10]針對常規(guī)層析方程系數(shù)計算求交運算量大的缺點，提出了一種提高運算速度和反演精度的投影面算法。

上述研究是利用完整穿過研究區(qū)域傾斜路徑上的水汽含量(slant water vapor,SWV)構(gòu)建層析模型的觀測方程，然而對于在層析區(qū)域側(cè)面穿出的射線，由于在組成觀測方程時并不能完全包含射線上的水汽信息，通常將該部分射線剔除。上述做法浪費了寶貴的GNSS觀測資源，也降低了觀測數(shù)據(jù)的利用率[11-12]。文獻[11]通過引入水汽密度比例因子的方法提高觀測數(shù)據(jù)的利用率，但該方法在計算水汽密度比例因子時需要探空數(shù)據(jù)提供的水汽密度參數(shù)作為支持。當(dāng)研究區(qū)域存在探空數(shù)據(jù)時，該方法有效可行；但若研究區(qū)域無探空數(shù)據(jù)時，則該方法會受到限制，因此，難以廣泛應(yīng)用。

針對傳統(tǒng)水汽層析方法在區(qū)域側(cè)面穿出射線無法利用的問題以及已有的側(cè)面射線利用方法中存在的缺陷，本文提出了利用ECMWF產(chǎn)品改善觀測數(shù)據(jù)利用率的方法，通過引入比例因子，使得在研究區(qū)域側(cè)面穿出的所有射線都能參與層析觀測方程的建立，大大提高了射線利用率和網(wǎng)格覆蓋率。此外，利用的ECMWF格網(wǎng)數(shù)據(jù)是全球分布的，能夠適用于全球任何區(qū)域。

1 利用ECMWF改善射線利用率的三維水汽層析算法

1.1 層析原理

GNSS衛(wèi)星信號穿過低空大氣層時，受到的斜路徑濕延遲(slant wet delay,SWD)可以表達成如下形式

(1)

利用式(2)即可得到斜路徑上的SWV

SWV=Π·SWD

(2)

(3)

式中，ρv表示水汽密度，單位為g/m3。根據(jù)地基GNSS對流層層析原理，將層析區(qū)域離散化成若干個獨立的三維單元網(wǎng)格，則衛(wèi)星信號路徑上的水汽含量可離散化成如下[4]

SWV=∑ijk(aijk·xijk)

(4)

式中，aijk表示射線在(i,j,k)網(wǎng)格內(nèi)的截距；xijk表示(i,j,k)網(wǎng)格內(nèi)的待估水汽密度值。

由于缺乏足夠的觀測數(shù)據(jù)以及層析區(qū)域位置的特定性，組成層析觀測方程的衛(wèi)星信號并非空間最優(yōu)分布，導(dǎo)致觀測方程的設(shè)計矩陣病態(tài)，致使三維水汽反演結(jié)果精度較差[4]。為了解決該問題，通常需要加入一些約束信息[2,4,16]。在水平方向上，本文根據(jù)水汽在空間上連續(xù)變化的特點引入水平約束[16]；在垂直方向上根據(jù)水汽隨高度的變化，利用負(fù)指數(shù)函數(shù)建立相鄰網(wǎng)格內(nèi)水汽的函數(shù)關(guān)系[17-18]。因此，傳統(tǒng)層析方法的模型可表示為

(5)

式中，A、H和V分別表示觀測方程、水平約束方程和垂直約束方程的系數(shù)矩陣。

1.2 利用ECMWF改善三維水汽層析射線利用率的方法

針對傳統(tǒng)水汽層析方法無法利用在研究區(qū)域側(cè)面穿出射線的缺陷，提出了利用ECMWF ERA-interim格網(wǎng)產(chǎn)品改善射線利用率的三維水汽反演算法。該方法首先利用ECMWF ERA-interim提供的分辨率為0.125°×0.125°的格網(wǎng)氣象資料，計算層析區(qū)域每個網(wǎng)格內(nèi)的水汽密度初值，然后通過引入比例因子(表示一條射線在層析區(qū)域內(nèi)部分的水汽含量與該射線路徑上總水汽含量的比值)，得到在側(cè)面穿出射線和側(cè)面交點的高與比例因子之間的關(guān)系，進而通過該比例因子關(guān)系得到在層析區(qū)域側(cè)面穿出射線在層析區(qū)域內(nèi)部分的水汽含量，進一步參與觀測方程的建立。該方法的優(yōu)點是能夠利用所有在層析區(qū)域側(cè)面穿出的射線參與層析模型中觀測方程的建立，大大改善了射線利用率和網(wǎng)格覆蓋度。

下面給出該方法的具體實現(xiàn)過程，具體流程見圖1。

圖1 利用ECMWF改善三維水汽反演射線利用率流程Fig.1 Flow chart of improving the utilization rate of satellite ray for three-dimensional water vapor tomography using ECMWF

(1) 利用ECMWF ERA-interim格網(wǎng)數(shù)據(jù)計算并插值出層析區(qū)域每個網(wǎng)格內(nèi)的水汽密度初值。由水汽狀態(tài)方程可得出[19]

(6)

式中，e表示水汽分壓，單位為hPa；T表示溫度，單位為K。

(2) 在水平方向上縮小層析區(qū)域，使原來在層析區(qū)域頂部穿出的射線在側(cè)面穿出。如圖2所示。

圖2 改善三維水汽層析射線利用率原理Fig.2 Schematic diagram of improving the utilization rate of satellite ray for three-dimensional water vapor tomography

(3) 計算比例因子。如圖2所示，射線OP的水汽含量值SWVOE可表示為如下形式

(7)

(8)

(9)

依次縮小層析區(qū)域，直到接收機位于縮小的層析區(qū)域之外為止。通過該方法，可以得到射線OP上多個高度對應(yīng)的比例因子。

(4) 對某一測站重復(fù)步驟(1)—(3)，直至該測站上所有在層析區(qū)域頂部穿出的射線對應(yīng)不同高度上的比例因子全部求出。

(5) 對層析區(qū)域內(nèi)的所有測站重復(fù)步驟(1)—(4)。通過分析比例因子和射線與側(cè)面交點的高之間的關(guān)系(如圖3所示)，建立如下比例因子模型

λ=a0+a1·h+a2·h2+a3·exp(h)

(10)

式中，a0-a3表示比例因子模型的系數(shù)，可以通過最小二乘方法求得；h為射線與側(cè)面交點的高程。需要說明，上述模型表達式是根據(jù)試驗時段內(nèi)每天UTC 00:00,06:00,12:00和18:00的數(shù)據(jù)通過擬合得到，后面將會詳細(xì)介紹。

圖3 比例因子與高程之間的關(guān)系Fig.3 Relationship between the scale factor and height

(11)

進一步可以求出射線OQ在層析區(qū)域內(nèi)上的水汽含量SWVOR

(12)

式中，SWVOQ表示射線OQ上總的水汽含量值。

(7) 利用側(cè)面穿出射線構(gòu)建觀測方程，最終得到本文提出方法的層析模型

(13)

式中，As表示側(cè)面穿出射線在研究區(qū)域網(wǎng)格內(nèi)的截距組成的系數(shù)矩陣；ys表示側(cè)面穿出射線在層析區(qū)域內(nèi)的水汽含量組成的列向量。

2 層析策略

2.1 試驗介紹

選取香港衛(wèi)星定位參考站網(wǎng)(satellite positioning reference station network,SatRef)中12個測站的數(shù)據(jù)為例進行分析，時間為2013年5月4日至5月21日共18 d的數(shù)據(jù)，各測站的地理位置如圖4所示。選取該時段數(shù)據(jù)進行試驗是因為該時段內(nèi)香港對應(yīng)著晴天、多云、陣雨、中雨和大雨等天氣情況，能夠很好地反映不同天氣下水汽的變化，具有一定的代表性。表1給出了試驗時段內(nèi)對應(yīng)的天氣狀況。

圖4 香港SatRef中GPS測站和無線電探空站地理位置分布Fig.4 Geographic distribution of GPS stations and Radiosonde station in Hong Kong SatRef

層析區(qū)域范圍：緯度22.19°N—22.54°N，經(jīng)度113.87°E—114.35°E，高度0～10 km；在緯向和經(jīng)向上的水平分辨率分別為0.05°和0.06°，垂直分辨率采用非等間距劃分方法[20]，從地面到層析頂部分別為0.5 km、0.5 km、0.7 km、0.7 km、0.7 km、0.9 km、0.9 km、0.9 km、1.1 km、1.1 km。因此，層析區(qū)域共有7×8×10個網(wǎng)格。此外，在研究區(qū)域內(nèi)還有一無線電探空站45004，如圖4中●所示。探空站每天在UTC00:00和12:00發(fā)射探空氣球，能夠觀測得到不同高度上的氣象數(shù)據(jù)，可作為本文水汽反演結(jié)果的檢驗。利用GAMIT軟件對采樣率為30 s的觀測數(shù)據(jù)進行處理，然后結(jié)合氣象數(shù)據(jù)得到各測站不同高度上的SWV[4,11]。

為了驗證本文方法，設(shè)計兩種方案進行層析試驗。方法1：利用傳統(tǒng)方法組成的層析模型反演水汽，如式(5)所示；方法2：利用本文方法組成的層析模型反演水汽，如式(13)所示。

表1 香港天氣統(tǒng)計

2.2 射線利用率及網(wǎng)格覆蓋率對比

本文首先對層析時段2013-05-04—2013-05-21共18 d的每天射線利用情況以及網(wǎng)格覆蓋率進行對比，如圖5所示。由圖可以明顯看出，利用本文方法(方法2)后，層析時段內(nèi)每天的射線利用率有了很大的提高，此外，網(wǎng)格覆蓋率也均有不同程度的提高。表2給出了18 d的平均射線利用率和有射線穿過的網(wǎng)格數(shù)的統(tǒng)計情況。通過計算可得，射線利用率平均提高了約54.55%，有射線穿過的網(wǎng)格數(shù)平均提高了16.43%，由原來的56.25%提高到72.68%。

圖5 試驗時段內(nèi)不同方法每天射線利用率及網(wǎng)格覆蓋率對比Fig.5 Comparison of utilization rate of satellite ray and coverage rate of voxels of each day derived from two method for the experimental period

表2 試驗時段內(nèi)射線利用情況和有射線穿過網(wǎng)格數(shù)統(tǒng)計結(jié)果

2.3 比例因子模型可靠性分析

由式(12)可以看出，層析區(qū)域側(cè)面穿出的射線在區(qū)域內(nèi)部分的水汽含量的精度依賴于比例因子模型的精度。因此，首先對試驗時段內(nèi)本文方法建立的比例因子模型的精度進行檢驗和可靠性分析。利用實測數(shù)據(jù)計算得到的比例因子與利用模型計算的比例因子進行對比，圖6分別給出了18 d比例因子模型的RMS和MAE。通過計算可得，層析時段內(nèi)比例因子模型的平均MAE為-0.007，這說明建立的模型系統(tǒng)誤差很??；該模型的平均RMS為0.054，通過統(tǒng)計，在研究區(qū)域側(cè)面穿出射線的SWV的均值約為110 mm，因此，由比例因子模型引起的誤差僅有6 mm左右，進一步驗證了建立的比例因子模型具有很高的精度。

圖6 試驗時段內(nèi)比例因子模型每天的RMS和MAE分布情況Fig.6 The RMS and MAE of scale factor model for the experimental period

3 試驗分析

3.1 內(nèi)符合精度檢驗

層析結(jié)果的質(zhì)量是檢驗建立層析模型精度的關(guān)鍵。因此，本文首先對上述兩種層析方法進行內(nèi)符合精度驗證。圖7給出了兩種層析方法計算得到的SWV殘差隨高度角變化的對比情況。由圖可以看出，兩種方法的SWV殘差隨高度角總體呈現(xiàn)遞減的趨勢，且方法2的SWV殘差波動比方法1的要小，這說明本文方法(方法2)其內(nèi)符合精度優(yōu)于傳統(tǒng)方法。通過計算，發(fā)現(xiàn)兩種方法的標(biāo)準(zhǔn)差(STD)和MAE分別為11.4/5.8 mm和10.6/4.4 mm，相對于傳統(tǒng)方法，本文方法的內(nèi)符合精度提高了約7.5%。

圖7 層析時段內(nèi)兩種層析方法計算得到的SWV殘差隨高度角變化對比圖Fig.7 The SWV residuals change with elevation derived from two tomographic method for the experimental period

3.2 外符合精度檢驗

前已述及，探空數(shù)據(jù)能夠提供精確的水汽密度廓線信息。此外，ECMWF ERA-interim也能夠提供高精度的分層氣象參數(shù)，因此，本文將這兩種不同來源數(shù)據(jù)計算的結(jié)果作為標(biāo)準(zhǔn)對水汽反演結(jié)果進行檢驗。首先利用不同水汽反演方法得到探空站所在位置上UTC00:00和12:00的水汽密度信息，然后分別與探空數(shù)據(jù)與ECMWF計算結(jié)果進行對比。圖8和9分別給出了試驗時段內(nèi)每天不同方法層析結(jié)果與探空數(shù)據(jù)和ECMWF對比的平均RMS和MAE。由兩圖可以看出，無論是與探空數(shù)據(jù)還是ECMWF對比，本文方法(方法2)得到的水汽質(zhì)量均優(yōu)于傳統(tǒng)方法(方法1)。表3給出了試驗時段內(nèi)不同方法層析結(jié)果與探空數(shù)據(jù)和ECMWF數(shù)據(jù)對比信息。由表可得，與探空數(shù)據(jù)對比，方法2的RMS和MAE分別為1.2 g/m3和0.3 g/m3，優(yōu)于方法1的1.7 g/m3和0.5 g/m3，其RMS的精度提高了29.4%；與ECMWF數(shù)據(jù)對比，方法2的RMS和MAE分別為2.1 g/m3和1.4 g/m3，優(yōu)于方法1的2.4 g/m3和1.6 g/m3，其RMS的精度提高了12.5%。

此外，對兩個特殊歷元(5-04UTC 00:00和5-21UTC 00:00)上的水汽廓線進行對比，選取該兩個歷元是因為它們對應(yīng)著試驗時段內(nèi)最小和最大的RMS值。圖10給出了兩個歷元上不同方法反演得到的水汽廓線與探空數(shù)據(jù)和ECMWF數(shù)據(jù)對比圖。由圖可以看出，兩種層析模型反演得到的水汽廓線與探空數(shù)據(jù)均有很好的一致性，ECMWF數(shù)據(jù)計算得到的水汽廓線精度稍差。通過計算，與探空數(shù)據(jù)、ECMWF數(shù)據(jù)對比，兩個歷元上方法1和方法2反演水汽密度的RMS分別為1.6/0.5 g/m3和1.9/1.5 g/m3。

圖8 試驗時段內(nèi)不同方法層析結(jié)果與探空數(shù)據(jù)對比每天的平均RMSFig.8 Comparison of average RMS of each day of water vapor density between Radiosonde，ECMWF and different tomographic result for the experimental period

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

圖9 試驗時段內(nèi)每天不同方法層析結(jié)果與探空數(shù)據(jù)對比的平均MAEFig.9 Comparison of average MAE of each day of water vapor density between Radiosonde，ECMWF and different tomographic result for the experimental period

圖10 兩個歷元不同方法反演得到的水汽廓線與探空數(shù)據(jù)和ECMWF數(shù)據(jù)對比Fig.10 Comparison of water vapor profile between Radiosonde，ECMWF and different tomographic result for two specific epochs

為了進一步對比不同方法反演的水汽密度與高程的關(guān)系，對層析時段內(nèi)不同高程上的水汽進行對比。圖11和圖12分別給出了兩種層析方法水汽結(jié)果在不同高度上與探空數(shù)據(jù)和ECMWF數(shù)據(jù)對比的RMS及相對誤差(relative error)。由兩圖可以看出，無論與探空數(shù)據(jù)還是ECMWF數(shù)據(jù)對比，在絕大多數(shù)高度上，方法2反演得到的水汽密度廓線的RMS和相對誤差均優(yōu)于方法1，這進一步證明了本文方法的優(yōu)越性。綜上對比可得，本文提出的三維水汽反演方法在香港區(qū)域不同天氣情況下均適用，因此，可以推斷該方法具有較好的實用性。

圖11 兩種層析方法水汽結(jié)果在不同高度上與探空數(shù)據(jù)對比的RMS和相對誤差Fig.11 RMS and Relative error changes with height derived from two tomographic result when compared with Radiosonde

圖12 兩種層析方法水汽結(jié)果在不同高度上與ECMWF數(shù)據(jù)對比的RMS和相對誤差Fig.12 RMS and Relative error changes with height derived from two tomographic result when compared with ECMWF

最后，圖13給出了2013年5月12日至14日對應(yīng)多云、晴、陣雨和大雨天氣的三維水汽層析結(jié)果。由上圖可以看出，在12日多云天氣下，水汽密度在不同層上分布，尤其是底層均多于13日的晴天天氣。而在13日和14日的陣雨和大雨天氣下，該地區(qū)在各層上的水汽密度值明顯高于多云和晴天天氣下的水汽分布，并且在14日大雨天氣下各層上的水汽密度值均大于13日的水汽反演結(jié)果，在給出的4 d中最大。綜合上述4 d不同天氣下不同層上的水汽密度分析可以看出，三維水汽密度值大小可依次排序為大雨>陣雨>陰天>晴天，這也與不同天氣下實際的水汽分布相符合。

圖13 2013年5月12日至14日水汽層析結(jié)果的三維分布圖Fig.13 Three-distribution of water vapor tomographic result during the May 12 to 14，2013

此外，本試驗選取的實驗時段(共18 d)對應(yīng)著晴天、多云、陣雨、中雨和大雨等多種天氣情況，在上述不同氣象條件下，水汽的垂直結(jié)構(gòu)分布有時并不符合負(fù)指數(shù)分布規(guī)律，但通過與探空數(shù)據(jù)對比以及對水汽的三維分布進行分析，發(fā)現(xiàn)本文方法仍然適用，均能取得可靠的反演結(jié)果。這是因為本文算法只是利用負(fù)指數(shù)分布函數(shù)對垂直方向上的網(wǎng)格進行約束，并不能決定最終的水汽反演結(jié)果。

4 結(jié) 論

本文提出了基于ECMWF產(chǎn)品改善三維水汽反演射線利用率的方法，通過引入比例因子，計算得到在層析區(qū)域側(cè)面穿出射線在層析區(qū)域內(nèi)部分的水汽含量，然后將其作為觀測方程輸入值參與水汽層析模型的建立。

利用香港SatRef中2013-05-04—2013-05-21共18 d的實測GPS和氣象數(shù)據(jù)對提出的方法進行驗證，通過對建立的比例因子模型進行分析，發(fā)現(xiàn)該模型計算SWV時的誤差僅為6 mm左右。本文方法能夠大大改善觀測數(shù)據(jù)的利用率以及射線穿過網(wǎng)格的覆蓋率，分別將探空數(shù)據(jù)和ECMWF數(shù)據(jù)計算結(jié)果作為標(biāo)準(zhǔn)對水汽反演結(jié)果精度及有效性進行驗證。試驗結(jié)果表明：本文方法能有效提高水汽反演結(jié)果的精度(分別提高了29.4%和12.5%)。此外，通過對不同高度上反演的水汽廓線進行對比，發(fā)現(xiàn)本文方法反演的水汽廓線的RMS和相對誤差在絕大多數(shù)層上均優(yōu)于傳統(tǒng)方法。

致謝：感謝IGRA和ECMWF提供的氣象數(shù)據(jù)，感謝香港SatRef提供的試驗數(shù)據(jù)。