韓 信，張寶忠，魏 征，李益農(nóng)，陳 鶴

（1.中國水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國家重點(diǎn)實驗室，北京 100038；2.中國農(nóng)業(yè)大學(xué) 水利與土木工程學(xué)院， 北京 100083）

1 研究背景

參考作物蒸騰量（ET0）是水資源與灌溉管理的重要數(shù)據(jù)，是水分平衡和灌溉調(diào)度的關(guān)鍵因素［1-3］，其計算和預(yù)測方法不僅成為農(nóng)田生態(tài)系統(tǒng)水分循環(huán)和水平衡研究的重要領(lǐng)域，同時在農(nóng)田灌溉制度優(yōu)化和水土資源配置等方面具有重要作用［4］。因此，為了更好地管理作物灌溉用水量以及提高作物水分利用效率，亟待對作物ET0進(jìn)行準(zhǔn)確預(yù)測。

1998年FAO-56推薦的Penman-Monteith（PM）模型是計算ET0的標(biāo)準(zhǔn)計算方法，能夠得到較精確的結(jié)果［1］。ET0與氣候變化等自然因素密切相關(guān)［5］。近年來涌現(xiàn)出諸多的ET0預(yù)測方法，如時間序列法［6-7］、灰色模型法［8］和經(jīng)驗公式法［9-10］。由于ET0與氣象因子之間存在復(fù)雜的非線性關(guān)系，常規(guī)的預(yù)測模型忽視其復(fù)雜的非線性關(guān)系造成預(yù)測精度不高，因此基于神經(jīng)網(wǎng)絡(luò)建立預(yù)測模型能夠考慮眾多因素對ET0的影響［11-12］，具有良好的應(yīng)用推廣價值。且徑向基人工神經(jīng)網(wǎng)絡(luò)（RBF）較Back Propagation（BP）神經(jīng)網(wǎng)絡(luò)具有更好的適用性［13-14］。Ozgur［15］用Levenberge Marquardt（LM）訓(xùn)練算法建立了基于RBF模型的ET0預(yù)測，并得到較好的結(jié)果，Sudheer等［16］利用有限的氣象資料建立了估算水稻日蒸散量的RBF模型，模型預(yù)測性能可靠。但是影響ET0的氣象因子本身也具有不確定性［17］，因此，近年來不確定性分析逐漸受到關(guān)注，但是考慮氣象因子不確定性的影響對ET0預(yù)測的研究較少。

本研究以華北平原農(nóng)田下墊面為研究對象，以研究區(qū)1981—2018年氣象因子和ET0（PM模型計算）資料序列為基礎(chǔ)，在貝葉斯理論框架下，根據(jù)馬爾科夫蒙特卡羅模擬與自適應(yīng)采樣算法結(jié)合的方法修正ET0影響因子，基于RBF模型構(gòu)建ET0不確定性預(yù)測模型。

2 預(yù)報方法簡介

2.1 研究區(qū)概況與數(shù)據(jù)來源 北京市大興區(qū)（39°26′—39°51′N，116°13′—116°43′E）位于華北平原北部永定河沖擊平原，總面積1031 km2，屬于溫帶半濕潤季風(fēng)氣候。季節(jié)變化明顯，夏季短而炎熱，冬季長而寒冷。多年平均氣溫12.1℃，平均最高氣溫14.8℃，最低氣溫6.7℃，平均相對濕度53%，年平均風(fēng)速1.84 m/s，年平均蒸散量1080 mm，多年年平均降雨量540 mm，7—9月降雨較多，降雨量占全年降雨總量的80%以上。北京市大興區(qū)下墊面主要為農(nóng)田，包括玉米/小麥、大豆，其中以玉米和小麥輪作為主，冬小麥整個生育期為260天左右（10月1日—次年6月30日），在正常年份冬小麥需補(bǔ)充灌溉，以保證作物對水分的需求。夏玉米生長期約90天左右（7月1日—9月30日），夏玉米全生育期不灌水［18］。

從中國氣象科學(xué)共享服務(wù)網(wǎng)（http：//cdc.cma.gov.cn）收集了北京市大興區(qū)1981—2018年的逐日歷史氣象數(shù)據(jù)和2013—2018年的逐日天氣預(yù)報數(shù)據(jù)，建模數(shù)據(jù)系列長度統(tǒng)一采用1981—2018年共38年的1—12月逐月數(shù)據(jù)。氣象數(shù)據(jù)包括最高溫度、最低溫度、平均溫度、平均相對濕度、最小相對濕度、日照時數(shù)、最大風(fēng)速、極大風(fēng)速、平均風(fēng)速、參考作物蒸散量。天氣預(yù)報數(shù)據(jù)及信息包括：最高溫度、最低溫度、相對濕度及天氣情況等，為了保持前后因子的一致性，最終選取最高溫度（Tmax）、最低溫度（Tmin）、相對濕度（RH）和平均風(fēng)速（Ws）建立ET0預(yù)測模型。天氣預(yù)報信息中Tmax、Tmin、RH可以直接利用，Ws則需要根據(jù)研究區(qū)天氣情況預(yù)報和氣象標(biāo)準(zhǔn)進(jìn)行解析［19］。根據(jù)風(fēng)速對地面的影響而引起的各種現(xiàn)象，按照風(fēng)速等級確定距離平面高處2 m的平均風(fēng)速（表1）。

表1 風(fēng)速等級

2.2 模型方法

2.2.1 參考作物蒸散量標(biāo)準(zhǔn)值 以從中國氣象科學(xué)共享服務(wù)網(wǎng)（http：//cdc.cma.gov.cn）獲取的研究區(qū)實測的逐日參考作物蒸散量（ET0）為標(biāo)準(zhǔn)值。

2.2.2 基于氣象因子不確定性的RBF模型構(gòu)建 利用徑向基人工神經(jīng)網(wǎng)絡(luò)模型（RBF-ANN）建立氣象因子與ET0的映射關(guān)系，并將其作為貝葉斯概率預(yù)報系統(tǒng)的似然函數(shù)，以ET0標(biāo)準(zhǔn)值為后驗信息，對ET0影響因子先驗信息進(jìn)行貝葉斯修正，利用基于自適應(yīng)算法的馬爾可夫鏈蒙特卡羅隨機(jī)模擬算法獲取各氣象因子的后驗密度。最后將修正后的因子結(jié)合RBF-ANN模型構(gòu)建考慮氣象因子不確定性的ET0預(yù)報模型（CU-RBF）。

（1）因子不確定性修正?；谪惾~斯概率預(yù)報系統(tǒng)對氣象因子的不確定性進(jìn)行修正。貝葉斯概率預(yù)報系統(tǒng)由Krzysztofowicz［20］提出，是一種通用的概率水文預(yù)報理論，可以與確定性水文模型耦合，考慮因子不確定性的ET0預(yù)測是在貝葉斯推斷的基礎(chǔ)上提出的，計算式為：

式中：h(y|x)為后驗分布；p(x|y)為x對y的條件概率；f(y)為y的邊緣密度，即先驗分布。

假設(shè)氣象數(shù)據(jù)X（x1為最高溫度、x2為最低溫度、x3為平均相對濕度、x4為平均風(fēng)速）的先驗密度為g(x)；ET0預(yù)測模型輸出參考作物蒸散量Sn已知時，X的似然函數(shù)為 f(S |X )。利用貝葉斯理論修正影響因子的不確定性，輸出Sn的期望值為：

式中：φ為X(X ={x1，x2， x3})的水文不確定性的定量表達(dá)；k為歸一化常數(shù)。

（2）修正后因子計算。AM算法是由Haario等［21］提出的一種有效的MCMC采樣器，其主要特點(diǎn)是隨機(jī)抽樣的參數(shù)空間初始協(xié)方差矩陣由先驗密度確定的多維正態(tài)分布，AM算法的主要步驟如下：

（a）初始化，i=0；

（b）計算轉(zhuǎn)移密度的協(xié)方差矩陣：

（3）CU-RBF模型構(gòu)建。徑向基函數(shù)網(wǎng)絡(luò)（RBF）由Broomhead提出［22］。該模型具有較快的收斂速度，在很多領(lǐng)域得到廣泛應(yīng)用［23-24］。RBF網(wǎng)絡(luò)模型主要包括輸入層、隱含層和輸出層，RBF模型的表達(dá)式為：

φ為可逆高斯函數(shù)：

式中：q為隱含層節(jié)點(diǎn)個數(shù)即樣本容量；wi為隱藏層各節(jié)點(diǎn)到輸出層F(X ) 的權(quán)因子；ci為第i個徑向基函數(shù)的中心值；δ為徑向基函數(shù)的變量，決定了該徑向基函數(shù)圍繞中心點(diǎn)的寬度。

模型的建立通過Matlab2017b的RBF工具箱來實現(xiàn)，具體函數(shù)的調(diào)用格式為：

式中：P、T分別為輸入和輸出變量；goal為均方誤差，取goal=0.001，避免過度擬合；spreads為RBF函數(shù)的寬度，通過隱含層神經(jīng)元的中心距離確定；MN為ANN學(xué)習(xí)過程中神經(jīng)元最大數(shù)目；DF為2次顯示之間添加神經(jīng)元的個數(shù)。

仿真函數(shù)調(diào)用為：

式中：Y為預(yù)測目標(biāo)值；X為輸入測試矩陣，net為式（9）中已經(jīng)學(xué)習(xí)好的網(wǎng)絡(luò)。

（4）訓(xùn)練集、驗證集及測試集。選擇1981—2018年1—12月444組修正后的歷史氣象月數(shù)據(jù)（Tmax、Tmin、RH和Ws）和以從中國氣象科學(xué)共享服務(wù)網(wǎng)（http：//cdc.cma.gov.cn）獲取的研究區(qū)逐日ET0累加得到的月ET0作為訓(xùn)練集；選擇2013—2015年36組修正后的歷史天氣預(yù)報月數(shù)據(jù)（根據(jù)天氣預(yù)報日值累加所得）和月ET0數(shù)據(jù)作為驗證集；選擇2016—2018年36組修正后的歷史天氣預(yù)報月數(shù)據(jù)（根據(jù)天氣預(yù)報日值累加所得）作為測試集。

2.3 模型評價 對于使用RBF神經(jīng)網(wǎng)絡(luò)模型預(yù)測值與PM計算標(biāo)準(zhǔn)值之間的相互比較，使用了5個評價指標(biāo)，包括決定系數(shù)（R2）、納什系數(shù)（NSE）、均方根誤差（RMSE）、平均絕對誤差（MAE）和相對誤差（RE）來評估模型的預(yù)測性能［25］，計算如下：

式中：xi為ET0預(yù)報值（利用修正后的氣象數(shù)據(jù)結(jié)合徑向基神經(jīng)網(wǎng)絡(luò)預(yù)測所得）；yi為ET0標(biāo)準(zhǔn)值（從中國氣象科學(xué)共享服務(wù)網(wǎng)（http：//cdc.cma.gov.cn）獲取的研究區(qū)逐日ET0）；y為標(biāo)準(zhǔn)值均值；i為預(yù)報樣本序列，i=1，2，…，x；n為預(yù)報值的樣本數(shù)。

3 結(jié)果與分析

3.1 數(shù)據(jù)集及先驗密度 圖1和表2為華北平原農(nóng)田下墊面1981—2018年氣象因子和ET0的數(shù)據(jù)系列，氣象因子包括最高溫度（Tmax）、 最低溫度（Tmin）、相對濕度（RH）、平均風(fēng)速（Ws）和ET0。Tmax、Tmin、RH、Ws和ET0多年平均值分別為15.15℃、6.74℃、44.20%、1.89 m/s和1063 mm，多年偏差系數(shù)分別為-0.36、-0.34、0.65、-0.36和-0.68，多年的氣象因子和ET0都在均值周圍有微小的波動，標(biāo)準(zhǔn)偏差分別為0.49℃、0.53℃、3%、0.16 m/s和0.1 mm。

表2 氣象參數(shù)和ET0

圖1 研究區(qū)1981—2018年氣象參數(shù)及ET0資料系列

圖2為氣象因子的正態(tài)概率分布。由圖2可知，各氣象因子均近似服從正態(tài)分布，通過SPSS軟件分析所知各因子的正態(tài)性假設(shè)均通過5%的顯著性水平。ET0影響因子的先驗密度和標(biāo)準(zhǔn)方差見表3，由表3可以看出，各月Ws的分布差異性較大，變差系數(shù)（CV均值為0.13）比Tmax（CV均值為0.11）、Tmin（CV均值為0.06）和RH（CV均值為0.07）的大，表明Ws相對其他氣象因子的不確定性較強(qiáng)。

表3 各月ET0影響因子的先驗密度

圖2 各氣象因子的正態(tài)概率圖

3.2 氣象因子不確定性識別 由于本研究所選氣象因子為4個變量，即Tmax、Tmin、RH、Ws，因此C0矩陣對角線元素取為各變量的先驗方差。設(shè)定AM-MCMC算法的初始抽樣次數(shù)為2000，運(yùn)行次數(shù)38次，每次采樣10 000個，這樣共抽取有效模擬樣本數(shù)為304 000個，用于每月氣象因子后驗分布的統(tǒng)計分析。氣象因子在迭代2000次后已經(jīng)收斂，所抽取的樣本已具有其總體的統(tǒng)計特征。

分別以3月（作物生長早期）、7月（作物生長中期）和10月（作物生長后期）典型月為例，所抽取各影響因子的后驗密度直方圖見圖3—圖5。Tmax、Tmin、RH、Ws的后驗密度分別為：最高溫度x1～N（10.46，0.982），最低氣溫x2～N（1.49，0.562），相對濕度x3～N（34.62，0.022），平均風(fēng)速x4～N（2.27，0.132）；最高溫度x1～N（25.73，0.522），最低氣溫x2～N（18.54，0.352），相對濕度x3～N（58.97，0.012），平均風(fēng)速x4～N（1.65，0.062）；最高溫度x1～N（16.03，0.682），最低氣溫x2～N（7.03，0.412），相對濕度x3～N（49.69，0.012），平均風(fēng)速x4～N（1.60，0.072）。對于Ws先驗密度的修正較Tmax、Tmin和RH大，在相同的顯著性水平條件下，各影響因子的置信區(qū)間不同，其中Ws比Tmax、Tmin和RH的范圍大，說明Ws具有較強(qiáng)的不確定性，主要是因為Ws受到外界的干擾較為敏感。因此，后驗信息對Ws的修正程度較其他3個因子大。Pereira等［4］的研究結(jié)果也表明為了避免因Ws的不確定性帶來的作物耗水預(yù)測不確定性，在進(jìn)行耗水預(yù)測時將不考慮Ws，而僅考慮T和RH等因子。由此可見，考慮Ws進(jìn)行ET0預(yù)測時，對其不確定進(jìn)行修正是必要的。

圖3 3月份頻率直方圖和密度曲線

圖4 7月份頻率直方圖和密度曲線

圖5 10月份頻率直方圖和密度曲線

3.3 參考作物蒸散量預(yù)測

3.3.1 單氣象因子不確定性修正對ET0預(yù)測影響 由3.2節(jié)可知，在相同的顯著性水平條件下，各氣象因子的均值置信區(qū)間范圍不同，其中Ws比Tmax、Tmin和RH的范圍大，說明Ws表現(xiàn)出較強(qiáng)的不確定性。圖6和表4為預(yù)測期單獨(dú)考慮修正某一因子（Tmax、Ws和RH）后RBF預(yù)測ET0的效果。

由圖6和表4可知，當(dāng)單獨(dú)對Tmax修正后，基于RBF模型的ET0的預(yù)測值與標(biāo)準(zhǔn)值的變化趨勢基本一致，且效果較未修正前好，模型評價參數(shù)效果均有所提高，R2、NSE分別增加0.04、0.01，RMSE、MAE和RE分別減少0.72 mm、0.02 mm和0.22%，當(dāng)分別單獨(dú)考慮Ws和RH修正后，基于RBF模型的ET0的預(yù)測值與標(biāo)準(zhǔn)值的變化趨勢基本一致，效果均較未修正前好，模型評價參數(shù)效果均有所提高。R2、NSE分別增加0.06、0.01和0.02、0.01，RMSE、MAE和RE分別減少6.08 mm、0.17 mm、1.71%和0.21 mm、0.01 mm、0.14%。由圖6對比可知，單獨(dú)對Ws修正，相較于分別單獨(dú)修正Tmax和RH時的預(yù)測效果較好，可能是因為風(fēng)速與氣象因子以外的農(nóng)田下墊面、觀測高度等有關(guān)，修正后的Ws能夠削弱外界條件對其的不確定性影響。而單獨(dú)對RH進(jìn)行修正時預(yù)測效果較差，R2、NSE分別相差0.05、0.02，RMSE、MAE和RE分別相差5.87mm、0.16mm、1.57%，可能是因為RH與降雨有直接關(guān)系，該研究結(jié)果與Mohmmad、刑貞相等研究結(jié)果一致［26-27］，在降雨條件下，各個氣象因子的不確定性增加，未對其他氣象因子進(jìn)行修正必然會造成預(yù)測的誤差。

圖6 最高溫度、相對濕度和平均風(fēng)速分別修正前后ET0預(yù)測效果

表4 分別單獨(dú)修正氣象因子不確定性修正前后ET0預(yù)測效果

3.3.2 考慮氣象因子不確定性修正的ET0預(yù)測 以AM-MCMC算法隨機(jī)模擬每月8000組ET0影響因子相應(yīng)的RBF預(yù)報值為對應(yīng)月份ET0的預(yù)測集合，將均值作為該月的ET0預(yù)測值，并得出指定概率（85%）的置信區(qū)間，實現(xiàn)了考慮了RBF模型輸入不確定性對ET0預(yù)報不確定性的影響（CU-RBF）。其中2013—2015年為驗證期，2016—2018年為測試期，驗證期和測試期的各月ET0計算與實測過程以及RBF的確定性和不確定性預(yù)報過程分別見圖7和圖8。

圖7 驗證期和測試期RBF、CU-RBF與標(biāo)準(zhǔn)值相關(guān)性分析

圖8 驗證期和測試期RBF、CU-RBF預(yù)測值與標(biāo)準(zhǔn)值對比

由圖7（a）、圖8（a）和表5可以看出，驗證期（2013—2015）和測試期（2016—2018）在未考慮氣象因子不確定性下建立的ET0預(yù)測模型（RBF）預(yù)測值與基于PM模型計算的標(biāo)準(zhǔn)值的相關(guān)性R2分別為0.87和0.89，擬合回歸系數(shù)分別為0.83和0.81，NSE分別為0.88、0.93。RBF確定性模型預(yù)測得到的ET0與PM模型計算的標(biāo)準(zhǔn)值趨勢基本一致，月均值分別為91.28 mm和88.98 mm，均方根RMSE和MAE分別為50.91 mm、46.27 mm和1.41 mm、1.29 mm；均值RE分別為6.99%、3.34%。由圖7（b）、圖8（b）和表5可以看出，驗證期和測試期在考慮氣象因子不確定性下建立的ET0概率預(yù)測模型（CURBF）預(yù)測值與基于PM模型計算的標(biāo)準(zhǔn)值的相關(guān)性的R2分別為0.94和0.97，擬合回歸系數(shù)分別為1.18和0.99，NSE分別為0.90、0.99。CU-RBF確定性模型預(yù)測得到的ET0與PM模型計算的標(biāo)準(zhǔn)值趨勢基本一致，相較于單獨(dú)考慮Tmax、RH、Ws的預(yù)測模型效果較好，NSE增加0.08；月均值分別為92.18 mm和91.05 mm，RMSE和MAE分別為47.28 mm、38.43 mm和1.31 mm、1.07 mm；均值RE分別為1.33%、-0.02%。綜上分析可知，無論在驗證期和測試期，CU-RBF模型預(yù)測的ET0與標(biāo)準(zhǔn)值的R2、NSE均有所提高，RMSE、MAE和RE均有所下降。R2、NSE分別提高0.07、0.08和0.09、0.06，RMSE、MAE和RE分別下降3.36 mm、7.84 mm，0.1 mm、0.22 mm和5.36%、3.36%。且CU-RBF能給出預(yù)報值的85%置信區(qū)間，其中ET0標(biāo)準(zhǔn)值包含在指定的概率區(qū)間內(nèi)的比例較高。Sudheer、邢貞相和德佳碩等［16，27-28］研究同樣表明氣象因子的不確定性必然影響作物需水預(yù)測的精度，在進(jìn)行基于神經(jīng)網(wǎng)絡(luò)模型對參考作物蒸散量預(yù)測前對氣象因子的不確定性修正能夠較大程度的提升其預(yù)測精度。本研究驗證了氣象因子不確定性修正在參考作物蒸散量預(yù)測方面實際應(yīng)用的可行性，可為農(nóng)田下墊面的未來水分管理提供科學(xué)依據(jù)。

表5 CU-RBF和RBF逐月ET0預(yù)測精度評價

4 結(jié)論

本研究通過分析華北平原區(qū)農(nóng)田下墊面氣象因子和ET0的變化情況，結(jié)合貝葉斯理論，利用基于自適應(yīng)采樣算法的AM-MCMC求解ET0影響因子（最高溫度、最低溫度、相對濕度、平均風(fēng)速）的后驗概率密度，應(yīng)用RBF-ANN建立ET0與氣象因子的函數(shù)關(guān)系，進(jìn)行ET0預(yù)報，得到以下結(jié)論：

（1）各氣象因子均近似服從正態(tài)分布，且Ws的變差系數(shù)比Tmax、Tmin和RH的大，相對其他主要?dú)庀笠蜃拥牟淮_定性較強(qiáng)。

（2）典型月份（3月、7月和 10月）的后驗密度分別為：Tmax～N（10.46，0.982），Tmin～N（1.49，0.562），RH～N（34.62，0.022），Ws～N（2.27，0.132）；Tmax～N（25.73，0.522），Tmin～N（18.54，0.352），RH～N（58.97， 0.012）， Ws～N（1.65， 0.062）； Tmax～N（16.03， 0.682）， Tmin～N（7.03， 0.412）， RH～N（49.69，0.012），Ws～N（1.60，0.072）。

（3）單獨(dú)修正Ws不確定性的CU-RBFWs模型比分別單獨(dú)修正Tmax和RH的CU-RBF預(yù)測模型好，R2、NSE分別為0.95、0.92，RMSE、MAE和RE分別為40.19 mm，1.12 mm和1.63%。

（4）在驗證期和測試期，考慮氣象因子不確定性構(gòu)建的CU-RBF模型的預(yù)報精度均有很大程度的提高。R2、NSE分別提高0.07、0.08和0.09、0.06，RMSE、MAE和RE分別下降3.36 mm、7.84 mm，0.1 mm、0.22 mm和5.36%、3.36%，且CU-RBF能得到85%置信區(qū)間的預(yù)報值。