楊 林，鄭南山，高淑照，劉 晨，李 楊，高士健

（1.西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院，成都 611756；2.中國(guó)礦業(yè)大學(xué) 環(huán)境與測(cè)繪學(xué)院，江蘇 徐州 221116；3.武漢大學(xué) 測(cè)繪學(xué)院，武漢 430072）

0 引言

水汽是地球大氣的組成成分之一，占大氣總體積的4%左右[1]，具有分布不均勻、變化頻率快的特點(diǎn)。它在三態(tài)變化過(guò)程中吸收或釋放大量熱量，是造成暴雨、洪澇等災(zāi)害的主要原因[2]；因此如何準(zhǔn)確預(yù)報(bào)水汽分布對(duì)于天氣預(yù)報(bào)具有重要的實(shí)用價(jià)值。

目前獲取大氣水汽的方法較多，如無(wú)線(xiàn)電探空技術(shù)、水汽輻射計(jì)技術(shù)等，但這些技術(shù)存在分辨率低、使用成本高的問(wèn)題。近年來(lái)隨著全球衛(wèi)星導(dǎo)航系統(tǒng)（global navigation satellite system，GNSS）的發(fā)展，可用GNSS 信號(hào)的天頂濕延遲（zenith wet delay，ZWD）分量及水汽轉(zhuǎn)換系數(shù)（Π）計(jì)算臺(tái)站上空的可降水量（precipitable water vapor，PWV），即地基GNSS 大氣水汽的反演[3]。利用該方法獲取水汽具有時(shí)空分辨率高、精度高、全天候的優(yōu)點(diǎn)[4]。高精度的ZWD 可通過(guò)GAMIT 軟件獲得。Π的計(jì)算涉及到大氣加權(quán)平均溫度Tm參數(shù)，因此能否準(zhǔn)確求解Tm值將直接決定PWV的反演質(zhì)量。文獻(xiàn)[1]建立了Tm與臺(tái)站溫度Ts的一元回歸方程。文獻(xiàn)[4]驗(yàn)證了全球定位系統(tǒng)（global positioning system，GPS）水汽反演在我國(guó)的可行性。文獻(xiàn)[5]進(jìn)行了大范圍的GPS 實(shí)驗(yàn)。1996 年，文獻(xiàn)[6]引入遠(yuǎn)距離基線(xiàn)來(lái)消除系統(tǒng)誤差，提高反演精度。文獻(xiàn)[7]證明地理位置和季節(jié)對(duì)Tm與Ts的相關(guān)性有重要影響。文獻(xiàn)[8]建立了適用于武漢地區(qū)的本地化加權(quán)平均溫度模型。文獻(xiàn)[9]證實(shí)了Tm與Ts存在非線(xiàn)性關(guān)系。文獻(xiàn)[10]證明本地化模型可更好地獲得高精度遙感水汽，但多因子模型相對(duì)單因子模型改進(jìn)效果不明顯。文獻(xiàn)[11]為了解決高海拔地貝維斯（Bevis）模型適用性較低的問(wèn)題，構(gòu)建了全球性Tm模型。文獻(xiàn)[12]構(gòu)建了貴州整體、局地及季節(jié)性Tm模型，證明其模型精度比Bevis 模型精度高。文獻(xiàn)[13]構(gòu)建了顧及精細(xì)季節(jié)變化的垂直遞減率函數(shù)模型，在此基礎(chǔ)上建立中國(guó)區(qū)域的Tm格網(wǎng)新模型。

在眾多的本地化Tm模型中，其建模方法多為基于最小二乘的線(xiàn)性回歸。文獻(xiàn)[9]提出Tm與Ts存在非線(xiàn)性關(guān)系?？紤]到反向傳播（back propagation，BP）神經(jīng)網(wǎng)絡(luò)算法對(duì)非線(xiàn)性關(guān)系具有強(qiáng)大的學(xué)習(xí)能力[14]，本文探討使用BP 神經(jīng)網(wǎng)絡(luò)算法建立一種能有效提升該地區(qū)PWV 反演質(zhì)量的Tm模型，并對(duì)模型的可行性進(jìn)行驗(yàn)證。

1 模型方法

在計(jì)算可降水量PWV 之前，需要先計(jì)算總延遲量即天頂對(duì)流層延遲（zenith tropospheric delay，ZTD），該值可在GNSS 數(shù)據(jù)處理的過(guò)程中獲得。ZTD 可分為干延遲（zenith hydrostatic delay，ZHD）和濕延遲ZWD。使用薩斯塔莫伊寧（Saastamoinen）模型可分離出高精度的ZHD[15]，其計(jì)算公式為

式中：ZHD 表示為ZHD；H為測(cè)站大地高，單位為mm；θ為測(cè)站緯度，單位為（°）；Ps為測(cè)站表面大氣壓，單位為hPa。該模型的誤差很小，精度一般可達(dá)毫米級(jí)。得到ZHD 后，即可計(jì)算ZWD 為

式中：ZWD 表示為ZWD；ZTD 表示為ZTD。然后將得到的ZWD與轉(zhuǎn)換系數(shù)Π相乘即可得到PWV，其計(jì)算公式為

式中：PWV 表示為PWV；Π為大氣水汽含量與濕延遲的無(wú)量綱比例因子，其計(jì)算公式為

式中：ρw為液態(tài)水密度，ρw=1×103kg/m3；Tm為大氣加權(quán)平均溫度，單位為K；e為沿垂直路徑上的水汽分壓，單位為hPa；T為大氣絕對(duì)溫度，單位為K；z0為積分起始高度，單位為m；Rv為水汽的氣體常數(shù)，Rv=461.495 J/kg；=（22.13±2.20）K/hPa；k3為大氣物理參數(shù)，k3=（3.739±0.012×105）K2/hPa。e、T為時(shí)空變化量，所以要想求得Π就必須先確定Tm值。式（4）雖能獲得高精度的加權(quán)平均溫度，但每天2 次的大氣探空探測(cè)不能提供連續(xù)的氣溫、氣壓等數(shù)據(jù)，無(wú)法完成積分運(yùn)算，所以常使用其離散化形式，即

式中：ei為第i層平均水汽壓，單位為hPa；Ti為第i層平均絕對(duì)值溫度值，單位為K；hi為第i層重力高程，單位為m；Td為露點(diǎn)溫度，單位為℃；rh為相對(duì)濕度，是一個(gè)百分量；ew為飽和水氣壓，是一個(gè)過(guò)渡參數(shù)，單位為hPa。用離散化形式計(jì)算的Tm值精度較高，可作為實(shí)驗(yàn)真值。

Bevis 等[1]在美國(guó)通過(guò)多次探空數(shù)據(jù)發(fā)現(xiàn)，地表溫度與大氣加權(quán)平均溫度之間存在緊密聯(lián)系，并根據(jù)探空資料建立線(xiàn)性回歸方程為

式中：Ts為地面溫度，單位為K；Tm為加權(quán)平均溫度，單位為K。由于式（6）的探空數(shù)據(jù)來(lái)源于中緯度的美國(guó)地區(qū)，直接用于中國(guó)則因存在地域差異性，所以很多地區(qū)未必適用；所以近年來(lái)越來(lái)越多的學(xué)者基于該方法使用本地?cái)?shù)據(jù)進(jìn)行建模。

2 加權(quán)平均溫度模型的建立

2.1 數(shù)據(jù)來(lái)源

本次實(shí)驗(yàn)的數(shù)據(jù)類(lèi)型包括2 類(lèi)：①香港國(guó)王公園探空站（King’s Park 45004）數(shù)據(jù)，每日的世界協(xié)調(diào)時(shí)（coordinated universal time，UTC）時(shí)間00:00 及12:00 時(shí)發(fā)射探空氣球，非極端天氣情況下每日產(chǎn)生2 組數(shù)據(jù)。選取2016 年、2017 年的探空數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。經(jīng)數(shù)據(jù)預(yù)處理后，2016 年得到653 個(gè)樣本數(shù)據(jù)，其中前70%的樣本數(shù)據(jù)用于模型的擬合，后30%用于模型的驗(yàn)證。2017 年得到的711 個(gè)樣本數(shù)據(jù)全部用于模型測(cè)試。②香港地區(qū)連續(xù)運(yùn)行參考站（continuously operating reference stations，CORS）HKSL的數(shù)據(jù)，選取00:00:00 UTC及12:00:00 UTC 時(shí)段用于計(jì)算對(duì)流層延遲分量。為減少對(duì)流層延遲相關(guān)性，引入3 個(gè)國(guó)內(nèi)的國(guó)際GNSS 服務(wù)組織（International GNSS Service，IGS），即北京房山站（BJFS）、上海佘山站（SHAO）和拉薩站（LHAZ），使用GAMIT 軟件一起進(jìn)行基線(xiàn)處理[16]，從而獲取HKSL 站的對(duì)流層延遲數(shù)據(jù)。IGS站及探空站的詳細(xì)信息如表1 所示，各站分布及站點(diǎn)間基線(xiàn)網(wǎng)情況如圖1 所示。

2.2 香港地區(qū)本地化模型的建立

文獻(xiàn)[9]中已介紹加權(quán)平均溫度與各氣象元素的線(xiàn)性關(guān)系不夠明確，而神經(jīng)網(wǎng)絡(luò)算法對(duì)于非線(xiàn)性映射關(guān)系具有強(qiáng)大的學(xué)習(xí)能力。BP 神經(jīng)網(wǎng)絡(luò)是一種導(dǎo)師制的學(xué)習(xí)網(wǎng)絡(luò)，它的工作原理是將一組訓(xùn)練集送入網(wǎng)絡(luò)，通過(guò)調(diào)整連接權(quán)，使網(wǎng)絡(luò)的實(shí)際輸出更加符合期望輸出。神經(jīng)網(wǎng)絡(luò)理論早已經(jīng)證明[17]：采用西格莫伊德（Sigmoid）轉(zhuǎn)換函數(shù)時(shí)，具有1 個(gè)隱層的BP 網(wǎng)絡(luò)可實(shí)現(xiàn)對(duì)任意函數(shù)的任意逼近。使用BP 神經(jīng)網(wǎng)絡(luò)模型進(jìn)行Tm模型的構(gòu)建，選取與Tm具有一定映射關(guān)系的4 個(gè)氣象參數(shù)作為輸入因子。這4 個(gè)氣象參數(shù)分別為露點(diǎn)溫度Td、測(cè)站壓強(qiáng)Ps、測(cè)站溫度Ts、水汽壓e。將離散積分計(jì)算的Tm值視為實(shí)驗(yàn)真值。

模型各階段收斂效果如圖2 所示。由圖可知：訓(xùn)練模型的均方誤差可達(dá)1×10-2數(shù)量級(jí)，可見(jiàn)模型擬合精度較高；通過(guò)對(duì)模型進(jìn)行驗(yàn)證，迭代12 次即可得到最佳收斂效果，可見(jiàn)模型收斂速度較快；最后對(duì)模型進(jìn)行測(cè)試，測(cè)試所得到的均方誤差同樣可以達(dá)到訓(xùn)練模型的數(shù)量級(jí)，可見(jiàn)訓(xùn)練模型的預(yù)測(cè)性能較好。

圖2 模型收斂效果

考慮到神經(jīng)網(wǎng)絡(luò)擬合難以得到直接的數(shù)學(xué)表達(dá)式，其訓(xùn)練后的關(guān)鍵參數(shù)仍保存在算法的隱含層中。為了方便對(duì)訓(xùn)練模型進(jìn)行精度評(píng)估，將模型的輸出值與真值進(jìn)行線(xiàn)性分析。圖3 展示了訓(xùn)練階段和測(cè)試階段模型歸一化輸出值與歸一化真值之間的線(xiàn)性關(guān)系。其中橫坐標(biāo)表示歸一化后輸入模型的Tm真值，即模型的學(xué)習(xí)目標(biāo)（目標(biāo)值），用t表示；縱坐標(biāo)表示歸一化后的Tm模型預(yù)測(cè)值，用p表示。根據(jù)相關(guān)系數(shù)的等級(jí)劃分，認(rèn)為0.7

圖3 模型擬合及測(cè)試情況

式中t為歸一化后的模型學(xué)習(xí)目標(biāo)值，即為T(mén)m的真值。在模型訓(xùn)練完成后，使用該模型對(duì)2017 年每日的UTC 時(shí)間00:00 和12:00 時(shí)的Tm值進(jìn)行測(cè)試，其預(yù)測(cè)效果與真值的線(xiàn)性關(guān)系如圖3（b）所示。由圖可知，測(cè)試階段真值與預(yù)測(cè)值之間的相關(guān)系數(shù)R=0.859。該值略低于訓(xùn)練樣本的相關(guān)系數(shù)，但其值仍在相關(guān)系數(shù)強(qiáng)相關(guān)性范圍內(nèi)，歸一化的Tm預(yù)測(cè)值與歸一化的Tm真值仍然呈現(xiàn)強(qiáng)相關(guān)性?？梢?jiàn)BP 神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)效果較好，也證明了BP 神經(jīng)網(wǎng)絡(luò)算法預(yù)測(cè)Tm值的可行性。測(cè)試階段模型歸一化輸出值與歸一化真值的線(xiàn)性表達(dá)式為

使用Bevis 模型對(duì)香港地區(qū)2017 年（UTC0Z和UTC12Z）的Tm值進(jìn)行計(jì)算。結(jié)合本地化模型的Tm預(yù)測(cè)值分別同Tm真值進(jìn)行比較（如圖4 所示）。由圖4（a）可知，Bevis 模型的預(yù)測(cè)值均小于真值，出現(xiàn)系統(tǒng)性的偏差；而圖4（b）本地化模型的預(yù)測(cè)值與真值基本一致。其原因是以上2 種模型的建模數(shù)據(jù)存在地區(qū)性的差異，可見(jiàn)Bevis 經(jīng)典模型不適用于香港地區(qū)的Tm預(yù)測(cè)，而本地化模型卻能很好地滿(mǎn)足本地區(qū)需求。

圖4 Bevis 模型和本地化模型的Tm 預(yù)測(cè)值與真值的比較

通過(guò)對(duì)2 種模型預(yù)測(cè)值的殘差進(jìn)行統(tǒng)計(jì)（如圖5 所示）：由圖5（a）可知，本地化模型的殘差主要分布在—2～2 K 之間，預(yù)測(cè)值的殘差在0 K 左右波動(dòng)，殘差分布較為均勻；圖5（b）中Bevis 模型的殘差主要分布在0～4 K 之間，整體偏大，可見(jiàn)本地化模型的預(yù)測(cè)性能更加穩(wěn)定。

圖5 不同模型的殘差分布

均方根誤差（root mean square error，RMSE）能反映預(yù)測(cè)值與真值之間的偏差，RMSE 值越小，預(yù)測(cè)精度越高。RMSE的計(jì)算公式為

式中：Xi為模型預(yù)測(cè)值；為真實(shí)值；N為樣本數(shù)。由表2 可知，本地化模型和Bevis 模型的RMSE分別為2.057 和3.236 K，本地化模型精度較Bevis模型提升了約36.4%，可見(jiàn)本地化模型所預(yù)測(cè)的加權(quán)平均溫度精度更高。

表2 不同模型計(jì)算Tm的誤差

3 本地化模型下PWV的反演分析

為了驗(yàn)證模型在實(shí)際應(yīng)用中的預(yù)測(cè)性能，首先使用GAMIT 軟件計(jì)算出香港地區(qū)CORS 站的ZTD，然后通過(guò)式（2）計(jì)算出ZWD，最后通過(guò)式（3）計(jì)算2 種模型下的PWV 預(yù)測(cè)值。將預(yù)測(cè)的PWV值與PWV 真值數(shù)據(jù)（來(lái)源于探空數(shù)據(jù)）進(jìn)行比較。圖6 為本地化模型和Bevis 模型對(duì)PWV 值的預(yù)測(cè)情況，圖中橫坐標(biāo)為樣本數(shù)，縱坐標(biāo)為PWV 值。綜合分析圖6 可得，2 種模型反演的大氣水汽含量與探空真值的趨勢(shì)基本一致。

圖6 不同模型反演的PWV 值與探空真值

依據(jù)式（9）計(jì)算2 種模型PWV 預(yù)測(cè)值的均方根誤差，結(jié)果如表3 所示：本地化模型和Bevis模型反演的PWV 與探空真值的均方根誤差分別為2.5 和3.0 mm，可見(jiàn)本地化模型的預(yù)測(cè)效果較Bevis模型提高了16.7%；2 種模型預(yù)測(cè)PWV的平均偏差分別為1.7 和2.4 mm，本地化模型的預(yù)測(cè)偏差較Bevis 模型降低了29.1%，可見(jiàn)本地化模型明顯提高了GNSS 大氣水汽反演精度。該實(shí)驗(yàn)結(jié)論與使用相對(duì)PWV 誤差評(píng)價(jià)Tm模型精度的結(jié)論一致，即本地化模型較Bevis 模型能更好地預(yù)測(cè)香港地區(qū)的PWV，該模型可作為香港地區(qū)預(yù)測(cè)高精度大氣水汽反演的可靠性模型。

表3 HKSL 站上不同Tm 模型反演PWV的誤差

4 結(jié)束語(yǔ)

本文利用香港地區(qū)2016 年的探空資料，使用BP 神經(jīng)網(wǎng)絡(luò)算法建立本地化的Tm模型，并對(duì)模型的精度進(jìn)行了分析，得到以下結(jié)論：①使用本地化模型對(duì)同地區(qū)2017 年的Tm值進(jìn)行預(yù)測(cè)，通過(guò)與Tm真值進(jìn)行比較，得到本地化模型的均方根誤差為2.057 K，較Bevis 模型降低約36%，可見(jiàn)本地化的Tm模型對(duì)香港地區(qū)的加權(quán)平均溫度預(yù)測(cè)性能優(yōu)于 Bevis 模型；②提取香港地區(qū)2017 年1 月的PWV 真值數(shù)據(jù)，與本地化模型的PWV 預(yù)測(cè)值進(jìn)行比較，得到其平均偏差和均方根誤差分別為1.7 和2.5 mm，平均偏差較Bevis模型降低了0.7 mm，均方根誤差較Bevis 模型降低了0.5 mm，可見(jiàn)本地化模型較Bevis 模型對(duì)香港地區(qū)大氣水汽反演有明顯的改善。但值得說(shuō)明的是，本文利用BP 神經(jīng)網(wǎng)絡(luò)建立本地化Tm模型僅基于本地?cái)?shù)據(jù)構(gòu)建，需要進(jìn)一步探究適用性更廣的模型。