全球氣候變化模式和多基因遺傳規劃在逐日參考作物騰發量預報中的應用

姚順秋，閆曉惠

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

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

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

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

1 研究方法

1.1 氣象數據

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

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

1.2 Hargreaves參考作物騰發量模型

Hargreaves 模型可表示為[12]:

(1)

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

Rs可通過下式計算[13]：

(2)

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

Ra的計算公式為[14]:

(3)

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

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

(4)

(5)

式中，J—日序號。

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

ωs=arccos(-tanφtanδ)

(6)

1.3 傳統與多基因遺傳規劃

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

1.4 誤差分析指標

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

(7)

(8)

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

2 研究結果與分析

2.1 基于GCM-MGGP法的ET0預報

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

圖1 MGGP模型的訓練演化過程

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

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

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

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

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

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

圖3 最優MGGP模型的樹狀結構圖

圖4 ET0的實際值與GCM-MGGP預報值

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

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

3 結語

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