基于MRI-CGCM3模式和遺傳編程人工智能算法的逐日參考作物騰發(fā)量預(yù)報(bào)方法

姚順秋，閆曉惠

(1.大連市莊河水利建筑勘測設(shè)計(jì)院，遼寧 大連 116400；2.加拿大渥太華大學(xué)，安大略 渥太華 K1N6N5)

作物騰發(fā)量是指作物葉面蒸騰量與顆間蒸發(fā)量之和，它的估算是作物需水量計(jì)算的基礎(chǔ)，也因此是精準(zhǔn)灌溉與排水設(shè)計(jì)的重要環(huán)節(jié)[1- 2]。作物騰發(fā)量可以表示為作物系數(shù)與參考作物騰發(fā)量(ET0)的乘積，而各類作物的作物系數(shù)基本為恒定值，因此，預(yù)報(bào)ET0是預(yù)報(bào)作物騰發(fā)量的關(guān)鍵、具有重要的現(xiàn)實(shí)意義[3- 5]。ET0的預(yù)報(bào)大體可以分為直接法和間接法兩種。直接法主要是根據(jù)歷史數(shù)據(jù)序列并采用基于數(shù)理統(tǒng)計(jì)原理的方法來對未來的ET0進(jìn)行預(yù)測。間接法主要是通過氣象預(yù)報(bào)數(shù)據(jù)并采用基于物理原理的ET0公式來對其進(jìn)行預(yù)測。應(yīng)用直接法可以對未來的演變趨勢進(jìn)行較為便捷迅速的判斷，但是，在全球氣候顯著變化的條件下，歷史數(shù)據(jù)的統(tǒng)計(jì)規(guī)律已難以完全切合未來的演變趨勢，因此直接法在逐日ET0的預(yù)測中已較難取得理想的精確度。

隨著近幾年氣象預(yù)測與模擬技術(shù)的發(fā)展，采用間接法預(yù)報(bào)逐日ET0正逐漸流行。例如閆曉惠等人[6]提出了基于全球氣候變化模式和Penman-Monteith公式的長期逐日ET0預(yù)報(bào)方法，并應(yīng)用該方法對加拿大渥太華、溫哥華、溫尼伯等城市的ET0進(jìn)行計(jì)算，結(jié)果證明基于該方法的預(yù)測值滿足精度要求。該方法的一個(gè)限制條件是它一般必須要求降尺度操作，操作復(fù)雜，且需要詳細(xì)的氣溫、日照、風(fēng)速和濕度歷史數(shù)據(jù)，在資料缺乏的地區(qū)較難推廣使用。近年來，人工智能與機(jī)器學(xué)習(xí)技術(shù)得到了飛速的發(fā)展，例如，閆曉惠等人[7- 11]采用將遺傳編程和多基因遺傳編程等方法應(yīng)用在了不同的水利問題中，取得非常好的效果。該方法主要具有3項(xiàng)重要優(yōu)勢：①可以降低人為操作的復(fù)雜性；②可以規(guī)避人為假設(shè)的不合理性；③可以深度挖掘數(shù)據(jù)本身之間的隱藏關(guān)系，從而可以更為準(zhǔn)確地表達(dá)數(shù)據(jù)之間的聯(lián)系。因此，本文提出基于MRI-CGCM3模式和遺傳編程人工智能算法的逐日ET0預(yù)報(bào)方法(以下簡稱為“GCM-GP法”)，首先基于實(shí)測氣溫?cái)?shù)據(jù)、采用Hargreaves公式計(jì)算了大連市莊河地區(qū)2011年7月1日至2020年3月31日間的ET0。分別采用MRI-CGCM3模式的原始模擬數(shù)據(jù)(以下簡稱為“GCM法”)和GCM-GP法計(jì)算了該地區(qū)的ET0，并與實(shí)際數(shù)據(jù)進(jìn)行比較，為提升ET0預(yù)報(bào)的效率與精確度提供參考。

1 方法

1.1 ET0計(jì)算方法

大連莊河位于遼寧省東南部，黃海北岸，是遼南地區(qū)重要的水源地，近幾年水資源壓力不斷上升，優(yōu)化灌溉排水方案對于該地區(qū)的水資源利用與保護(hù)成為極為重要的措施，而ET0的準(zhǔn)確預(yù)報(bào)又是灌溉排水規(guī)劃的最重要環(huán)節(jié)之一，因此該地區(qū)亟需提高ET0的預(yù)測精度。獲取該莊河地區(qū)2011年7月1日至2020年3月31日的逐日最高與最低氣溫?cái)?shù)據(jù)，數(shù)據(jù)共3197組。該地區(qū)采用Hargreaves 公式[12]計(jì)算ET0，該公式可以表示為:

(1)

式中，T—日最高氣溫與最低氣溫的平均值，℃；Rs—太陽輻射，MJ/(m2·d)。

Rs可通過下式計(jì)算[13]：

(2)

式中，KRs—經(jīng)驗(yàn)系數(shù)，對于內(nèi)陸地區(qū)其值一般設(shè)定為0.16，而對于沿海地區(qū)其值一般設(shè)定為019，Tmax、Tmin—日最高和最低氣溫；Ra—地外輻射，MJ/(m2·d)。

Ra的計(jì)算公式為[14]:

(3)

式中，GSC—太陽常數(shù)，為 0.0820;dr—日地相對距離;ωs—日落時(shí)角;φ—維度;δ—太陽偏磁角。

日地相對距離dr和太陽偏磁角δ的計(jì)算公式為：

(4)

(5)

式中，J—日序號。

日落時(shí)角ωs的計(jì)算公式為：

ωs=arccos[-tanφtanδ]

(6)

1.2 MRI-CGCM3模式

MRI-CGCM3是由氣象研究中心(Meteorological Research Institute: MRI)在2012年開發(fā)的一款全球氣候模式。它主要是在第5代耦合模式比較計(jì)劃(CMIP5)框架下開發(fā)，在MRI-CGCM2系列版本上進(jìn)行更新與改進(jìn)而成。該模式主要可提供氣溫(包括最高、最低和平均氣溫)、降雨、海平面氣壓、風(fēng)速和降雪5項(xiàng)氣象模擬數(shù)據(jù)。目前，在各項(xiàng)氣象因子的預(yù)測中，氣溫的預(yù)測精確度一般較高，而Hargreaves公式中只需要逐日氣溫?cái)?shù)據(jù)，因此，選取其中的逐日最高和最低氣溫。采用MATLAB程序提取莊河地區(qū)所在位置(39.6808°N，122.9673°E)的數(shù)據(jù)。

1.3 遺傳編程

遺傳編程也稱基因編程，它是一種利用進(jìn)化算法的機(jī)器學(xué)習(xí)技術(shù)，它主要是受達(dá)爾文進(jìn)化論的啟發(fā)，是借鑒生物演化過程而產(chǎn)生的一種可以構(gòu)造算法的算法。它首先隨機(jī)生成模型庫，之后對模型進(jìn)行性能評價(jià)，再通過基因繁殖(Reproduction)、基因突變(Mutation)和基因交叉(Crossover)等運(yùn)算來對模型進(jìn)行演化，最終得到滿意的模型。本文采用該算法來建立ET0模型并進(jìn)行ET0的計(jì)算，不同于傳統(tǒng)的模型，在該方法中，輸入量為MRI-CGCM3模式的氣溫?cái)?shù)據(jù)，輸出量為實(shí)際的ET0，因此采用該算法得到的最終模型可以直接表示MRI-CGCM3模式氣溫?cái)?shù)據(jù)與實(shí)際ET0之間的關(guān)系，省卻了降尺度與地區(qū)修正等操作過程，因此，模型建立后的使用非常簡便。

1.4 性能評價(jià)

為深入量化預(yù)報(bào)方法的精確度，計(jì)算模擬值與實(shí)測值之間的均方根誤差(RMSE)和決定系數(shù)值(R2)，其公式如下。

均方根差：

(7)

決定系數(shù)：

(8)

式中，xs—實(shí)測值；xm—模擬值。

2 結(jié)果與討論

2.1 逐日氣溫實(shí)測與模擬值

提取并整理研究時(shí)間范圍內(nèi)逐日最高與最低氣溫的實(shí)測與MRI-CGCM3模擬值(以下簡稱“模擬值”)，結(jié)果呈現(xiàn)于圖2中，其中圓點(diǎn)表示實(shí)測值，實(shí)線代表模擬值。由圖可知，模擬逐日氣溫隨時(shí)間的變化趨勢與實(shí)測值基本保持一致，各散點(diǎn)的分布較接近于實(shí)測值。實(shí)測逐日最高氣溫?cái)?shù)據(jù)序列中的最低值與最高值分別為-19°和35℃，而模擬值數(shù)據(jù)序列中的最低值與最高值分別為-20℃和33℃，兩組數(shù)據(jù)相差較小。實(shí)測逐日最低氣溫?cái)?shù)據(jù)序列中的最低值與最高值分別為-23°和26℃，而模擬值數(shù)據(jù)序列中的最低值與最高值同樣也分別為-23℃和26℃，與實(shí)測數(shù)據(jù)完全一致。因此，MRI-CGCM3模擬值較為合理，可以用于后續(xù)的計(jì)算與分析中。

圖1 逐日最高與最低氣溫實(shí)測值與模擬值

2.2 基于GCM法的ET0預(yù)報(bào)結(jié)果

圖2比較了基于實(shí)測數(shù)據(jù)的ET0計(jì)算結(jié)果(以下簡稱“實(shí)際值”)和基于GCM法的計(jì)算結(jié)果(以下簡稱“GCM模擬值”)。大多數(shù)的各數(shù)據(jù)點(diǎn)均分布在等值線附近，說明GCM模擬值與實(shí)際值的吻合度較好，但部分散點(diǎn)超過了20%的范圍。多數(shù)超過20%范圍的散點(diǎn)位于等值線的下側(cè)，即采用當(dāng)前的方法可能低估實(shí)際的ET0值。因此，該方法的精確度有待提高。GCM模擬結(jié)果的RMSE值為1.099毫米/天、R2值為0.76。

圖2 逐日ET0的實(shí)測值與GCM模擬值

2.3 基于GCM-GP的ET0預(yù)報(bào)結(jié)果

采用MATLAB程序從3197組數(shù)據(jù)中隨機(jī)選取80%用于模型訓(xùn)練，剩余20%的數(shù)據(jù)用于模型驗(yàn)證。從MRI-CGCM3模式中得到逐日最高與最低氣溫?cái)?shù)據(jù)并定義為輸入量，選用實(shí)際的ET0值為輸出量，采用該算法將得到一個(gè)新的ET0模型，之后可以采用該模型和輸入量計(jì)算輸出量。該算法首先隨機(jī)生成一個(gè)能夠擬合特定數(shù)據(jù)集的函數(shù)庫(此處設(shè)定為500個(gè)函數(shù))，因此在初代時(shí)其誤差較大(平均RMSE值約為1.25)。隨著進(jìn)化代數(shù)的增加，該算法通過基因繁殖、基因突變和基因交叉等運(yùn)算來對函數(shù)模型進(jìn)行演化與改進(jìn)，因此模型的誤差值不斷減小(圖3)。根據(jù)圖3可知，當(dāng)代數(shù)達(dá)到約50時(shí)，其誤差隨代數(shù)的變化已經(jīng)較小，即運(yùn)行更多的進(jìn)化代數(shù)將不再顯著提升模型性能。因此，將最大進(jìn)化代數(shù)設(shè)置為300，而最終得到的最佳函數(shù)模型可以簡化表示為：

y=0.00149 log(x1-1.0 absx2)2+0.00149 abs(10.0x1-5.0x2+sin(x1+2.0x2)+log(sin2.0x1)+28.5)+0.00149 abs(3.0x1-1.0x2+sin(x1+2.0x2)-14.2)-0.00298 sinx1+0.00149(x1+x2+21.2)(2.0x1-2.0x2+abs(3.0x1-2.0x2)+7.11)-0.00149x1(3.0x1-2.0x2-1.0 log(x1-2.0x2+7.03)2+abs(3.0x2-2.0x1+0.463)+7.11)-0.00149 abs(cos(x2-1.0x1+7.15)+log(x2/(x2+7.11))(8.0x1-7.0x2+sinx1-35.5)-0.108

式中，y—輸出量(ET0)；x1—輸入變量1(日最高氣溫)；x2—輸入變量2(日最低氣溫)。

圖3 模型訓(xùn)練演化過程

采用該模型計(jì)算ET0值(以下稱為GCM-GP模擬值)，并與實(shí)際值進(jìn)行對比，如圖4所示。

其中用于模型訓(xùn)練的數(shù)據(jù)序列的RMSE和R2值分別為0.37mm/d和0.93，而用于模型驗(yàn)證的數(shù)據(jù)序列的RMSE和R2值分別為0.37mm/d和0.94。兩者相差不大，說明模型性能穩(wěn)定。相對于GCM結(jié)果，其RMSE值降低了約66%，而R2值提高了約24%，因此可以認(rèn)為該方法可以使得ET0的精確度大幅度提升、具有非常好的推廣價(jià)值。

圖4 ET0的實(shí)際值(菱形)與GCM-GP模擬值(實(shí)線)

3 結(jié)語

采用莊河地區(qū)2011年7月1日至2020年3月31日逐日氣溫觀測數(shù)據(jù)與Hargreaves模型計(jì)算了該地區(qū)的逐日ET0。提出了基于MRI-CGCM3模式與遺傳編程人工智能算法(GCM-GP)的ET0估算方法，并分別用兩種方法對ET0進(jìn)行預(yù)報(bào)。

結(jié)果表明，兩種方法所得結(jié)果的變化趨勢基本與實(shí)測結(jié)果保持一致，但GCM-GP方法可以顯著提高ET0的模擬精確度(相對于GCM法，其RMSE值降低了約66%，而R2值提高了約24%)，具有非常好的推廣價(jià)值。本文在實(shí)際ET0的計(jì)算中采用的是Hargreaves模型，但是在資料比較充足的地區(qū)，Penman-Monteith公式的使用則更為普遍，因此，在下一步的研究中，可以采用相同的方法訓(xùn)練出對應(yīng)于Penman-Monteith公式的ET0預(yù)報(bào)模型，也可以采用該方法建立其它地區(qū)的ET0模型。