趙利江，楊海鵬，許超鈐，趙健赟

1. 長安大學(xué) 地質(zhì)工程與測繪學(xué)院，西安 710054；

2. 青海省基礎(chǔ)測繪院，西寧 810016；

3. 青海省地理空間信息技術(shù)和應(yīng)用重點(diǎn)實(shí)驗(yàn)室，西寧 810016；

4. 武漢大學(xué) 測繪學(xué)院，武漢 430079；

5. 青海大學(xué) 地質(zhì)工程系，西寧 810016

1 引 言

湟水流域是青海省政治、經(jīng)濟(jì)、文化和交通中心，流域內(nèi)人口數(shù)為312 萬，占全省人口數(shù)的60.2%；工農(nóng)業(yè)總產(chǎn)值為143 億元，占全省的54.2%。20 世紀(jì)90 年代以來隨著氣候增暖、降水減少、地表蒸發(fā)量的不斷增加，湟水流域氣候干旱化程度進(jìn)一步加快，異常天氣（雪災(zāi)、干旱、洪澇、冰雹等）發(fā)生次數(shù)明顯增加（戴升等，2006）。尤其是2022年8 月發(fā)生的持續(xù)降雨，引發(fā)了大通縣的山洪和泥石流，造成4 人死亡，27 人失聯(lián)。而水汽是降水發(fā)生的最直接因素，因此，研究湟水流域的水汽變化及其與強(qiáng)降水的相互關(guān)系，對指導(dǎo)防災(zāi)減災(zāi)工作具有十分重要的現(xiàn)實(shí)意義。

可降水量（precipitable water vapor，PWV）是從地面直到大氣頂界的單位面積大氣柱中所含水汽總量全部凝結(jié)并降落到地面可以產(chǎn)生的降水量，是反映大氣水汽含量主要技術(shù)指標(biāo)。近年來，隨著國家衛(wèi)星導(dǎo)航定位基準(zhǔn)站網(wǎng)和各省衛(wèi)星導(dǎo)航定位基準(zhǔn)站網(wǎng)的逐步建立，基于全球?qū)Ш叫l(wèi)星系統(tǒng)（global navigation satellite system，GNSS）反演PWV逐步成為研究熱點(diǎn)（馬進(jìn)全等，2019；張鵬等，2018；王洪棟，2018）。PWV精度受大氣加權(quán)平均溫度（Tm）影響。Tm可以使用探空站的大氣垂直氣象要素計(jì)算，但受限于站點(diǎn)數(shù)量較少且時間分辨率較低，難以滿足實(shí)際應(yīng)用需求，因此，傳統(tǒng)研究通常用Tm數(shù)據(jù)建立模型，此類模型主要分為回歸模型和經(jīng)驗(yàn)?zāi)Ｐ蛢煞N?；貧w模型通過分析Tm與地表氣象要素的相關(guān)關(guān)系采用回歸分析等方法建立，最早的是利用北美 13 個探空站兩年的探空記錄建立的Bevis 式模型（Bevis 等，1992）。但因其系數(shù)具有季節(jié)和地理的限制，在其他地區(qū)使用誤差較大。為提高本地模型精度，有研究基于我國探空數(shù)據(jù)建立了東部地區(qū)的Bevis 式模型（李建國等，1999；劉焱雄等，2000）；有研究利用Tm與地表氣溫、氣壓、露點(diǎn)溫度、水汽壓相關(guān)關(guān)系，建立了本地的多因子模型，結(jié)果表明增加因子數(shù)能夠提高Tm的估計(jì)精度（劉旭春等，2006；周國君和潘雄，2006；李國翠等，2008）。上述回歸模型計(jì)算Tm均需要地表氣象要素，難以滿足實(shí)時水汽的計(jì)算需要。因此，不需要?dú)庀髤?shù)的經(jīng)驗(yàn)?zāi)Ｐ完懤m(xù)出現(xiàn)，如GWMT-D（He 等，2017）、全球氣壓溫度3（global pressure and temperature 3，GPT3）（Landskron 和Bohm，2018）、GGTm（Huang 等，2019）等，其中，GPT3 模型具有開源、易操作、格網(wǎng)分辨率和精度高的特點(diǎn)，是目前使用最為廣泛的經(jīng)驗(yàn)?zāi)Ｐ汀PT3 模型出現(xiàn)之后，將積日（day of year，DOY）引入?yún)^(qū)域加權(quán)平均溫度模型，所得精度優(yōu)于 GPT3（莫智翔等，2021）。楊飛等（2022）優(yōu)化后的GPT3 模型中誤差達(dá)到3.52 K。但經(jīng)驗(yàn)?zāi)Ｐ途韧陀诨貧w模型（Ding，2020）。隨著大數(shù)據(jù)的發(fā)展，有研究利用神經(jīng)網(wǎng)絡(luò)基于氣溫、氣壓、水汽壓、高程、緯度、積日等建立了大氣加權(quán)平均溫度模型（謝劭峰等，2022）。相比傳統(tǒng)模型精度有一定的提升，卻未考慮大氣加權(quán)平均溫度的長期變化趨勢。

為進(jìn)一步研究提高大氣加權(quán)平均溫度模型的精度，本文利用歐洲中期天氣預(yù)報(bào)中心（European Centre for Medium-range Weather Forecasts ，ECMWF）的全球氣候第五代再分析（the fifth generation ECMWF atmospheric reanalysis，ERA5）數(shù)據(jù)集和大氣逐層數(shù)據(jù)，考慮大氣加權(quán)平均溫度的年際變化，基于多層感知器（multilayer perceptron，MLP）方法建立了湟水流域大氣加權(quán)平均溫度模型，并與已有的Bevis 式、雙因子、多因子、GPT3、改進(jìn)的GPT3 模型、謝劭峰等（2022）方法六種模型進(jìn)行了比較驗(yàn)證。結(jié)果表明，本文所建立的改進(jìn)的MLP 模型具有更高的精度。

2 數(shù)據(jù)來源

湟水流域只有西寧1 座探空站，因此僅使用探空數(shù)據(jù)建立的Tm模型雖然在局部有很高的精度，但無法滿足整個流域的計(jì)算需求。本文從ECMWF 收集了湟水流域2010~2020 年的大氣逐層數(shù)據(jù)、ERA5 數(shù)據(jù)集中1950 年以來的陸地地表氣溫?cái)?shù)據(jù)，主要包括大氣逐層氣壓、溫度、位勢高度、比濕和地表溫度、氣壓、露點(diǎn)溫度等參數(shù)。在國家氣象科學(xué)數(shù)據(jù)中心（http://data.cma.cn/）收集了11 座氣象站的地表氣溫、氣壓、露點(diǎn)溫度等數(shù)據(jù)。從美國懷俄明大學(xué)官網(wǎng)（http://www.weather.uwyo.edu/）下載了2010~2022 年西寧站探空數(shù)據(jù)，其中包括氣壓、高度、地面溫度、露點(diǎn)溫度、比濕、相對濕度等參數(shù)，用于檢驗(yàn)?zāi)Ｐ途?。另外，收集了流域?nèi)30 m 分辨率的數(shù)字高程模型，用于獲取格網(wǎng)平均高程。為了彌補(bǔ)氣象站數(shù)量的不足，使用陸地地表氣溫?cái)?shù)據(jù)補(bǔ)充了6 個格網(wǎng)點(diǎn)。氣象站及格網(wǎng)點(diǎn)分布，如圖1 所示。

圖1 站點(diǎn)分布Fig.1 Distribution map of stations

3 研究方法

3.1 Tm 計(jì)算原理

在利用GNSS 觀測數(shù)據(jù)進(jìn)行水汽反演過程中，首先要對GNSS 原始數(shù)據(jù)進(jìn)行處理以獲取天頂總延遲（zenith total delay，ZTD），并用獲取的估值減去由于大氣引起的天頂靜力學(xué)延遲（zenith hydrostatic delay，ZHD），從而得到與水汽相關(guān)的天頂濕延遲（zenith wet delay，ZWD）?？山邓颗c天頂濕延遲之間的關(guān)系如下：

式中，PWV 為可降水量，mm；Π為轉(zhuǎn)換因子，其與大氣加權(quán)平均溫度Tm函數(shù)關(guān)系為

式中，Rv為水汽氣體常數(shù)；k2、k3為常數(shù)。Rv=461.522 J/（kg·K），k2′=22.1±2.2，k3=3.739（±0.012）×105K2/hPa。

大氣加權(quán)平均溫度Tm與大氣水汽壓ei、大氣溫度Ti、層高dh的關(guān)系可使用ECMWF 提供的逐層大氣參數(shù)離散積分：

式中，ei、Ti分別為第i層大氣的平均水汽壓（hPa）和平均氣溫（K）；hi為第i層大氣厚度，m。hi由勢高作差得到，由于層間勢高差與高程差相差較小，因此，本文在數(shù)據(jù)處理時，直接使用勢高差代替高程差。由于ECMWF 提供的大氣逐層數(shù)據(jù)、陸地地表氣溫?cái)?shù)據(jù)中沒有水汽壓，需要使用比濕由式（4）或露點(diǎn)溫度由式（5）求得

3.2 多層感知器

MLP 是一種前向結(jié)構(gòu)的人工神經(jīng)網(wǎng)絡(luò)，由多個節(jié)點(diǎn)層所組成，其中，第一層為輸入層，最后一層為輸出層，中間部分為一個或多個隱含層。除輸入節(jié)點(diǎn)外，每個節(jié)點(diǎn)都是一個帶有非線性激活函數(shù)的多輸入單輸出神經(jīng)元。結(jié)構(gòu)如圖2 所示。

圖2 MLP 結(jié)構(gòu)示意Fig.2 Architecture of multi-layer perceptron network

3.3 精度評價指標(biāo)

使用平均偏差Sbia和均方根誤差Srm表示新模型與參考值間的離散程度：

式中，N為匹配樣本總量；yi為模型估計(jì)Tm值，K；zi為探空數(shù)據(jù)實(shí)測Tm值，K。

相關(guān)系數(shù)Corr 表示兩個樣本的相關(guān)關(guān)系：

式中，X為被檢驗(yàn)數(shù)據(jù)，包括地表溫度（K）、水汽壓（hPa）；Y為Tm樣本數(shù)據(jù)，K。

4 實(shí)驗(yàn)與分析

4.1 Tm 變化特征

由大氣逐層數(shù)據(jù)、陸地地表氣溫?cái)?shù)據(jù)中提取了湟水流域11 座氣象站、6 個格網(wǎng)點(diǎn)的136562 個樣本數(shù)據(jù)，計(jì)算了全部樣本數(shù)據(jù)的Tm值。圖3（a）描述的是2010~2020 年湟水流域的日均Tm隨時間的變化情況，湟水流域年均Tm約270 K，最高點(diǎn)約280 K（DOY 為210），最低點(diǎn)約250 K（DOY 為365），具有明顯的季節(jié)性變化。圖3（b）描述的是2010~2020 年湟水流域的年均Tm隨時間變化情況，除了季節(jié)性變化之外，該流域Tm年均值具有一定的年際變化：2010~2012 年下降；2012~2016 年回升；2016~2020 年下降。因此，本文將積日和累計(jì)日（cumulative day，CD）作為建模參數(shù)。

圖3 湟水流域2010~2020 年大氣加權(quán)平均溫度變化情況Fig.3 Changes of atmospheric weighted average temperature in Huangshui basin from 2010 to 2020

4.2 相關(guān)性分析

為明確Tm模型參數(shù)，本文計(jì)算了2010~2020 年湟水流域136562 個樣本數(shù)據(jù)的大氣加權(quán)平均溫度Tm，分析了Tm與地表氣溫Ts、水汽壓es、高程hs、緯度bs的相關(guān)性，其中，Ts與Tm相關(guān)系數(shù)0.934，es與Tm相關(guān)系數(shù)0.85，均呈正向強(qiáng)相關(guān)，因此，建模時應(yīng)考慮Ts、es。hs、bs對Tm有系統(tǒng)性的影響，如不考慮緯度和高程的變化則模型中會存在一定的系統(tǒng)誤差，Tm與Ts、es、hs、bs的相關(guān)關(guān)系，如圖4 所示。

圖4 大氣加權(quán)平均溫度與地表氣溫（a）、水汽壓（b）、高程（c）和緯度（d）的相關(guān)關(guān)系Fig.4 Correlation between atmospheric weighted average temperature with surface temperature (a), water vapor pressure (b),elevation (c), and latitude (d)

4.3 模型結(jié)果對比

4.3.1 典型模型對比

為了解不同Tm模型的精度水平，本文利用2010~2020 年湟水流域的ERA5 數(shù)據(jù)建立了三種有氣象要素的線性回歸分析模型（Bevis 式、雙因子、多因子）、一種經(jīng)驗(yàn)?zāi)Ｐ停℅PT3）、一種混合模型（改進(jìn)的GPT3）。并利用均方根誤差、平均誤差比較和評價了五種傳統(tǒng)模型估計(jì)Tm的精度，如表1 所示。

表1 傳統(tǒng)模型的參數(shù)及其內(nèi)符合精度Tab.1 Traditional model parameters and its internal coincidence accuracy K

（1）在Bevis 式模型的基礎(chǔ)上，增加水汽壓可使回歸模型均方根誤差減小17%，進(jìn)一步增加高程和緯度因子可降低回歸模型的平均誤差，但均方根誤差沒有明顯降低。

（2）經(jīng)驗(yàn)?zāi)Ｐ虶PT3 在湟水流域具有較好的適用性，在不需要地表氣象參數(shù)的條件下平均誤差可達(dá)0.078 K，標(biāo)準(zhǔn)差2.17 K。

（3）混合模型的均方根誤差優(yōu)于GPT3 和多因子模型，但平均誤差有明顯增大，需要使用探空數(shù)據(jù)進(jìn)一步驗(yàn)證。

4.3.2 改進(jìn)的MLP 模型

為了進(jìn)一步提高模型精度，借鑒謝劭峰等（2022）方法的多層感知器原理，本文引入?yún)?shù)——年際變化因子（CD），并結(jié)合地表氣溫Ts、水汽壓es、高程hs、緯度bs和DOY 建立了改進(jìn)的MLP 模型。

（1）選取與Tm具有相關(guān)關(guān)系的Ts、es、hs、bs、DOY、CD 作為協(xié)變量輸入到輸入層，Tm作為因變量。

（2）將70%數(shù)據(jù)作為訓(xùn)練集，30%作為驗(yàn)證集，反向優(yōu)化模型參數(shù)，利用2018 年西寧探空站實(shí)測Tm評估所建模型的性能。

（3）定義神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)，利用試湊法得到模型隱藏層的層數(shù)為1、節(jié)點(diǎn)數(shù)為4。隱藏層激活函數(shù)為雙曲正切函數(shù)，輸出層激活函數(shù)為恒等函數(shù)y=x，損失函數(shù)選用誤差平方和。

（4）在模型訓(xùn)練中，選取批次訓(xùn)練，優(yōu)化算法選擇標(biāo)度共軛梯度法。

研究中起點(diǎn)時間為2010 年1 月1 日。相比于謝劭峰等（2022）方法，本模型平均偏差減小50%，均方根誤差減小4%，如表2 所示。

表2 MLP 相關(guān)的模型參數(shù)及內(nèi)符合精度Tab.2 Model parameters and internal coincidence accuracy based on MLP K

5 精度評價

5.1 基于探空數(shù)據(jù)的模型精度分析

為客觀檢驗(yàn)?zāi)Ｐ途龋褂?010~2022 年西寧探空站實(shí)測Tm對Bevis 式、雙因子、多因子、GPT3、改進(jìn)的GPT3 模型、謝劭峰等（2022）方法與改進(jìn)的MLP 模型進(jìn)行了評價，結(jié)果如表3 所示?？紤]了CD因子的本模型均方根誤差相比Bevis 式、雙因子、多因子、GPT3、改進(jìn)的GPT3、謝劭峰等（2022）方法分別減小了32%、23%、15%、14%、7%、5%。此外，本模型顧及了地表溫度、水汽壓、高程、緯度、積日、累計(jì)日多種影響因素，模型表現(xiàn)得更穩(wěn)定，如圖5 所示。因此，本模型可以為青海湟水流域GNSS反演PWV 提供更高精度的大氣加權(quán)平均溫度值。

表3 研究所用模型外符合精度Tab.3 Out-of-model coincidence accuracy K

圖5 模型殘差散點(diǎn)分布Fig.5 Scatter plot of model residuals

5.2 Tm 對PWV 的影響估計(jì)

為量化Tm對PWV 的影響，利用式（1）、式（2）分別對П和Tm微分得式（9）、式（10）。Tm取湟水流域年均大氣加權(quán)平均溫度270 K，進(jìn)一步簡化得式（11）。如

式中，σП為轉(zhuǎn)換因子的均方根誤差；σPWV為PWV均方根誤差，mm。

式中，mTσ為Tm均方根誤差，K。

取湟水流域平均濕延遲150 mm。得到各模型計(jì)算PWV 精度情況，如表4 所示?？芍狟evis 式、雙因子、多因子、GPT3、改進(jìn)的GPT3、謝劭峰等（2022）方法、本模型的Tm誤差帶給PWV 的影響分別為0.336 mm、0.293 mm、0.267 mm、0.263 mm、0.245 mm、0.240 mm、0.227 mm。

表4 不同方法產(chǎn)生的PWV 誤差Tab.4 PWV error of different models

6 結(jié) 論

考慮了年際變化因素，本文利用青海湟水流域2010~2020 年大氣逐層數(shù)據(jù)和陸地地表氣溫?cái)?shù)據(jù)，建立了適用于青海湟水流域的改進(jìn)的MLP大氣加權(quán)平均溫度模型，并以2010~2022 年西寧探空站獲取的Tm為參考，與已有研究方法進(jìn)行了精度比較評價，結(jié)果表明，平均偏差和均方根誤差分別為–0.01 K、2.71 K，均方根誤差相比于Bevis 式、雙因子、多因子、GPT3、改進(jìn)的GPT3模型、謝劭峰等（2022）方法分別減小了32%、23%、15%、14%、7%、5%。這證明本文改進(jìn)的MLP 模型在青海湟水流域的精度要優(yōu)于已有研究方法，一定程度上能夠?yàn)殇宜饔虼髿庠u估研究提供方法參考。