全球氣候變化模式和多基因遺傳規(guī)劃在逐日參考作物騰發(fā)量預(yù)報中的應(yīng)用

姚順秋，閆曉惠

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

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

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

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

1 研究方法

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

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

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

(2)

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

Ra的計算公式為[14]:

(3)

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

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

(4)

(5)

式中，J—日序號。

日落時角ωs的計算公式為：

ωs=arccos(-tanφtanδ)

(6)

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

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

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

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

(7)

(8)

式中，xs—實測值；xm—模擬值。

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

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

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

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

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

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

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

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

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

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

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

圖4 ET0的實際值與GCM-MGGP預(yù)報值

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

圖6 逐日ET0的實際值與GCM-SPPG法預(yù)報值的對比散點圖

3 結(jié)語

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