姚順秋，閆曉惠

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

估算作物需水量是農(nóng)業(yè)水利工程規(guī)劃與設(shè)計(jì)中的關(guān)鍵環(huán)節(jié)[1- 2]，而作物需水量的估算主要是基于參考作物騰發(fā)量(ET0)的計(jì)算。因此，ET0的預(yù)報(bào)對(duì)水資源的精細(xì)化配置與管理工作具有非常重要的意義[3]。ET0的估算方法中較為簡(jiǎn)單的是直接估算法，即對(duì)歷史ET0資料進(jìn)行數(shù)據(jù)分析，并以此推測(cè)未來的ET0變化過程。近年來，全球氣候變化顯著，歷史數(shù)據(jù)規(guī)律已經(jīng)很難準(zhǔn)確描述未來的水文氣象變化趨勢(shì)，因此直接估算法難以滿足當(dāng)前的精準(zhǔn)化水資源管理要求。

間接估算法的原理是采用基于實(shí)際物理原理的ET0計(jì)算模型和氣象預(yù)報(bào)數(shù)據(jù)進(jìn)行ET0的計(jì)算，得益于氣象預(yù)測(cè)能力的不斷提高，間接估算法目前得到了大量的關(guān)注與研究。例如，劉夢(mèng)等[4]采用天氣預(yù)報(bào)數(shù)據(jù)對(duì)漳河灌區(qū)的參考作物騰發(fā)量進(jìn)行了預(yù)報(bào)研究，徐俊增等[5]也對(duì)基于天氣預(yù)報(bào)的參考作物騰發(fā)量預(yù)報(bào)方法進(jìn)行了比較。近期，閆曉惠等[6]采用Penman-Monteith模型、全球氣候變化模式和降尺度法對(duì)加拿大渥太華、溫哥華等6個(gè)城市的逐日ET0進(jìn)行了預(yù)報(bào)，取得良好的預(yù)報(bào)精確度。但是，該方法一方面對(duì)數(shù)據(jù)要求較高，另一方面也需要較為繁雜的降尺度操作，限定了其廣泛的應(yīng)用。隨著人工智能與機(jī)器學(xué)習(xí)技術(shù)的不斷進(jìn)步，采用人工智能算法來推演模型成為重要的技術(shù)手段。近年來，閆曉惠等[7- 9]成功將人工智能與機(jī)器學(xué)習(xí)技術(shù)引入到不同的水資源相關(guān)應(yīng)用中，且證明該項(xiàng)技術(shù)在水利領(lǐng)域中具有非常明顯的優(yōu)勢(shì)。首先，采用機(jī)器學(xué)習(xí)技術(shù)代替人為分析可以大幅度節(jié)省人力成本；其次，人工智能算法不需要提供預(yù)設(shè)的模型結(jié)構(gòu)，從而可以避免預(yù)設(shè)模型結(jié)構(gòu)的不合理性；此外，該方法也可以深度挖掘變量之間的相互作用關(guān)系，從而可提供更為精確的模型[10]。在各類機(jī)器學(xué)習(xí)算法中，遺傳編程在水利領(lǐng)域中具有廣闊的應(yīng)用前景，它的一個(gè)重要特點(diǎn)是可以提供一個(gè)明確的數(shù)學(xué)模型。多基因遺傳規(guī)劃(MGGP)是近期在傳統(tǒng)的遺傳編程算法上演變而來，相對(duì)于傳統(tǒng)遺傳編程方法，它具有精確度更高、復(fù)雜度更低等優(yōu)點(diǎn)，因此具有非常大的發(fā)展?jié)摿7- 9，11]。

但據(jù)筆者所知，目前，基于多基因遺傳編程和全球氣候變化模式的ET0估算方法幾乎沒有被報(bào)道過。因此，本文旨在對(duì)該方法進(jìn)行可行性研究與性能評(píng)測(cè)。首先采用大連市莊河地區(qū)2011年7月1日—2020年3月31日間的逐日最高與最低氣溫?cái)?shù)據(jù)、基于Hargreaves模型計(jì)算了該地區(qū)與時(shí)間段內(nèi)的逐日ET0值。分別采用本文新提出的方法(GCM-MGGP)、全球氣候變化模式的原始數(shù)據(jù)法(GCM法)和全球氣候變化模式與傳統(tǒng)的遺傳規(guī)劃方法(單基因遺傳規(guī)劃；GCM-SGGP)對(duì)該ET0序列進(jìn)行計(jì)算，并對(duì)各項(xiàng)結(jié)果進(jìn)行比較分析，為ET0的預(yù)報(bào)工作提供新的方法與經(jīng)驗(yàn)。

1 研究方法

1.1 氣象數(shù)據(jù)

實(shí)測(cè)氣象數(shù)據(jù)選取大連市莊河地區(qū)2011年7月1日—2020年3月31日共3197組數(shù)據(jù)的逐日最高與最低氣溫?cái)?shù)據(jù)。日最高氣溫?cái)?shù)據(jù)序列的最大值、最低值、平均值和中間值分別為35、-19、14.7、16°C；標(biāo)準(zhǔn)差、方差、峰度和偏度值分別為11.2、126.5、1.8、-0.3。日最低氣溫?cái)?shù)據(jù)序列的最大值、最低值、平均值、和中間值分別為26、-23、5.3、6°C；標(biāo)準(zhǔn)差、方差、峰度、和偏度值分別為11.8、138.4、1.8、-0.1。

全球氣候變化氣象預(yù)測(cè)數(shù)據(jù)主要提取自MRI-CGCM3模式。該模式是在CMIP5(第五代耦合模式比較計(jì)劃)框架下開發(fā)的全球氣候模式，主要包含氣溫、降雨、海平面氣壓、風(fēng)速和降雪5項(xiàng)氣象因子。本文提取其中對(duì)應(yīng)于莊河地區(qū)的模擬結(jié)果，其坐標(biāo)為北緯39.6808°、東經(jīng)122.9673°。在該數(shù)據(jù)中，日最高氣溫?cái)?shù)據(jù)序列的最大值、最低值、平均值和中間值分別為33.3、-20.1、11.0、11.7°C；標(biāo)準(zhǔn)差、方差、峰度、和偏度值分別為12.0、142.9、1.9、-0.3。日最低氣溫?cái)?shù)據(jù)序列的最大值、最低值、平均值和中間值分別為25.9、-23.6、5.3、5.8°C；標(biāo)準(zhǔn)差、方差、峰度和偏度值分別為12.0、144.7、1.9、-0.2。

1.2 Hargreaves參考作物騰發(fā)量模型

Hargreaves 模型可表示為[12]:

(1)

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

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

(2)

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

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

(3)

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

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

(4)

(5)

式中，J—日序號(hào)。

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

ωs=arccos(-tanφtanδ)

(6)

1.3 傳統(tǒng)與多基因遺傳規(guī)劃

遺傳規(guī)劃是基于達(dá)爾文進(jìn)化論和孟德爾遺傳變異理論思想、參考生物演進(jìn)過程而開發(fā)的一種可以構(gòu)造算法與模型的算法。該方法可以隨機(jī)產(chǎn)生模型種群，并對(duì)各模型進(jìn)行評(píng)測(cè)，若不滿足要求，則自動(dòng)采用基因繁殖、基因突變、和基因交叉等運(yùn)算來對(duì)種群進(jìn)行改進(jìn)，直至出現(xiàn)滿足要求的模型。在傳統(tǒng)的遺傳規(guī)劃算法中，一個(gè)模型染色體只包含有一個(gè)基因(因此也可稱為單基因遺傳規(guī)劃)，但在多基因遺傳規(guī)劃算法中，一個(gè)模型染色體可以有多個(gè)基因，從而可以使得其精確度更高或復(fù)雜度更低。在本項(xiàng)目的模型訓(xùn)練過程中，選取的輸入量分別為GCM最高氣溫和GCM最低氣溫，選取的輸出量為實(shí)際的ET0值。因此，所得模型可以直接建立GCM模擬氣溫?cái)?shù)據(jù)與實(shí)際ET0值之間的關(guān)系，而不需要進(jìn)行ET0模型計(jì)算、降尺度分析、和地區(qū)修正等操作，因此采用訓(xùn)練后的模型進(jìn)行ET0預(yù)報(bào)具有操作簡(jiǎn)便、易于上手的現(xiàn)實(shí)優(yōu)點(diǎn)。

1.4 誤差分析指標(biāo)

采用均方根誤差(RMSE)和決定系數(shù)值(R2)來量化預(yù)報(bào)值與實(shí)際值之間的誤差，其公式分別為：

(7)

(8)

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

2 研究結(jié)果與分析

2.1 基于GCM-MGGP法的ET0預(yù)報(bào)

將數(shù)據(jù)組隨機(jī)分配為兩個(gè)部分，分別為訓(xùn)練數(shù)據(jù)序列和驗(yàn)證數(shù)據(jù)序列。其中，訓(xùn)練數(shù)據(jù)占數(shù)據(jù)總量的80%，主要用于訓(xùn)練模型；剩余數(shù)據(jù)為驗(yàn)證數(shù)據(jù)，主要是當(dāng)做為未知數(shù)據(jù)以評(píng)測(cè)所得模型的預(yù)報(bào)性能。圖1呈現(xiàn)的為MGGP模型的訓(xùn)練演化過程。在第一代模型種群中，各模型為隨機(jī)產(chǎn)生，因此其誤差較大，RMSE值超過1mm/d，之后，MGGP算法利用演化運(yùn)算來對(duì)模型進(jìn)行改進(jìn)，只需10代左右，模型種群的誤差值變化幅度即已較小，說明運(yùn)行更多的進(jìn)化代數(shù)不再顯著提高種群的精確度。因此，將最終的總進(jìn)化代數(shù)設(shè)置為300即可滿足要求。

圖1 MGGP模型的訓(xùn)練演化過程

每代種群包含500個(gè)模型，圖2繪制的為最后一代種群中各模型的復(fù)雜度與性能指標(biāo)值。判斷一個(gè)模型的優(yōu)劣一般要考慮精確性和簡(jiǎn)易性兩個(gè)方面，而這兩個(gè)方面卻通常是矛盾的。本研究采用Pareto優(yōu)化法，并將位于Pareto優(yōu)化解曲線上的模型用圓圈表示。這些模型的特點(diǎn)是：在相同的精確度下，這些模型最為簡(jiǎn)易；而在同樣的復(fù)雜度下，這些模型的結(jié)果最為精確。最終，選擇這些模型中精確度最高的模型作為整個(gè)種群中的最優(yōu)解。最優(yōu)模型的樹狀結(jié)構(gòu)呈現(xiàn)于圖3中，其中x1代表日最高氣溫、x2代表日最低氣溫。

圖2 MGGP模型種群中各模型的復(fù)雜度與性能值

應(yīng)用該最優(yōu)MGGP模型計(jì)算研究區(qū)域的逐日ET0值，并稱之為“GCM-MGGP預(yù)報(bào)值”。圖4對(duì)比了ET0的實(shí)際值與GCM-MGGP預(yù)報(bào)值。由圖可知，GCM-MGGP數(shù)據(jù)與實(shí)際數(shù)據(jù)的變化規(guī)律基本保持一致。訓(xùn)練數(shù)據(jù)的RMSE值為0.365mm/d，R2值為0.936；驗(yàn)證數(shù)據(jù)的的RMSE值為0.364mm/d，R2值為0.938。RMSE值均較低，且R2值均較高，說明該方法可以提高滿意的預(yù)報(bào)精確度。同時(shí)，訓(xùn)練數(shù)據(jù)序列與驗(yàn)證數(shù)據(jù)序列的誤差水平非常接近，說明模型訓(xùn)練過程中的過擬合風(fēng)險(xiǎn)較低。

2.2 與GCM法和GCM-SGGP方法的比較

圖5呈現(xiàn)了研究區(qū)域內(nèi)逐日ET0的實(shí)際值與GCM法預(yù)報(bào)值的時(shí)間序列。由圖可知，GCM法大體上可以準(zhǔn)確地預(yù)報(bào)出數(shù)據(jù)的起伏變化過程，但是，較多的數(shù)據(jù)點(diǎn)誤差較大。而且，多數(shù)誤差較大的數(shù)據(jù)預(yù)報(bào)值低于實(shí)際值，不利于水資源管理的安全性。GCM法預(yù)報(bào)數(shù)據(jù)的總體RMSE值和R2值分別為1.099mm/d和0.76，因此，精確度較低。根據(jù)RMSE的關(guān)系可知，GCM-MGGP法相對(duì)于GCM法可以將誤差降低約67%，成效顯著。GCM-MGGP法可大幅度提高其預(yù)報(bào)精度的一個(gè)重要原因是它相當(dāng)于自動(dòng)增加了高精準(zhǔn)度的數(shù)據(jù)降尺度與地區(qū)修正操作并考慮了變量之間深度隱藏的作用關(guān)系。

圖6為逐日ET0的實(shí)際值與GCM-SPPG法預(yù)報(bào)值的對(duì)比散點(diǎn)圖。圖中實(shí)線為1∶1等值線。當(dāng)散點(diǎn)接近于等值線時(shí)，說明該數(shù)據(jù)點(diǎn)接近于實(shí)際值，否則說明偏差較大。圖6顯示，大部分散點(diǎn)均較接近于實(shí)際值，說明該方法的預(yù)報(bào)結(jié)果較接近實(shí)際值。較多的點(diǎn)位于等值線的左上方，說明該方法傾向于高估實(shí)際的ET0值。訓(xùn)練數(shù)據(jù)的RMSE值為0.373mm/d，R2值為0.934；驗(yàn)證數(shù)據(jù)的的RMSE值為0.368mm/d，R2值為0.936。RMSE值均高于GCM-MGGP法的結(jié)果，而R2值均低于GCM-MGGP法的結(jié)果，說明GCM-MGGP法相對(duì)于GCM-SGGP法更為精確。此外，采用Smits與Kotanchek[15]的Expressional-Complexity法衡量所得模型的復(fù)雜度，得到GCM-SGGP最優(yōu)解的復(fù)雜度為1797。如圖3所示，采用MGGP法所得的模型中各個(gè)基因?yàn)榈途S度線性或非線性項(xiàng)，而這些項(xiàng)是通過線性方向進(jìn)行組合，因此，模型的非線性維度隨基因數(shù)的增加而傾向于降低，本例中最優(yōu)模型的復(fù)雜度為766，遠(yuǎn)低于GCM-SGGP最優(yōu)解的復(fù)雜度。因此，可以認(rèn)為GCM-MGGP法相對(duì)于GCM-SGGP法不但可以提高預(yù)報(bào)精確度，也可以降低所得模型的復(fù)雜度。

圖3 最優(yōu)MGGP模型的樹狀結(jié)構(gòu)圖

圖4 ET0的實(shí)際值與GCM-MGGP預(yù)報(bào)值

圖5 逐日ET0的實(shí)際值與GCM法預(yù)報(bào)值的時(shí)間序列

圖6 逐日ET0的實(shí)際值與GCM-SPPG法預(yù)報(bào)值的對(duì)比散點(diǎn)圖

3 結(jié)語

分別采用GCM法、GCM-SGGP法、和本文新提出的GCM-MGGP法計(jì)算了大連市莊河地區(qū)2011年7月1日—2020年3月31日間的逐日ET0值。結(jié)果顯示，相對(duì)于GCM法，GCM-MGGP法可顯著提高預(yù)報(bào)精確度。相對(duì)于GCM-SGGP法，GCM-MGGP法預(yù)報(bào)精度也有所提高，而其Expressional-Complexity復(fù)雜度可大幅降低。因此，GCM-MGGP法預(yù)報(bào)精確、模型簡(jiǎn)易，具有較好的發(fā)展?jié)摿?。受限于觀測(cè)資料的不足，莊河地區(qū)實(shí)際ET0值的計(jì)算是基于Hargreaves模型，隨著水文氣象資料觀測(cè)能力的提升，未來可以采用更符合現(xiàn)實(shí)情況的Penman-Monteith模型進(jìn)行ET0計(jì)算，并采用類似方法訓(xùn)練出GCM-MGGP人工智能模型。此外，也可以將該方法在更多地區(qū)進(jìn)行應(yīng)用與驗(yàn)證，以進(jìn)一步評(píng)測(cè)其性能。